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CHAPTER  I.  INTRODUCTION  AND 
EXECUTIVE  SUMMARY 


Hft 


During  1980,  JAYCOR  has  continued  its  modeling  of  the  implosion 
dynamics  and  radiative  output  of  diode-imploded  annular  plasmas,  in 
close  coordination  with  work  ongoing  at  the  Naval  Research  Laboratory. 

This  report  treats  three  areas  of  advance  during  the  1980  effort: 

(I)  improvements  to  the  1-D  strongly-coupled  plasma  implosion  and  radia¬ 
tion  code  SPLAT  and  results  of  radiative  yield  studies  using  the  code; 

(II)  development  of  formalism  for  solving  the  field  penetration/skin- 
depth  problem  in  an  inhomogeneous,  time-varying  imploding  conductor  in  a 
plasma-loaded  diode;  (III)  circuit  equation  and  scaling  of  hard  radiation 
in  the  presence  of  fully  developed  sausage  instability  (beading)  of  the 
assembled  plasma.  In  addition,  a  short  section  (Chapter  V)  is  devoted  to 
work  in  progress:  high-accuracy  matrix  inversion  techniques  and  inter¬ 
polators  for  solving  the  generalized  Hertz  vector  equations  used  in  II 
above,  and  for  following  CRE  equations  and  diffusive  behavior  in  general; 
and,  beginning  plans  for  modifying  1-D  MHD  codes,  making  them  compatible 
with  the  field-diffusion  and  corona  programs  and  with  CRE  radiation 
packages . 

I.  The  code  Improvements  to  SPLAT  are: 

(A)  a  rudimentary  two-bin  spectral  resolution  of  the  radiative 
output  (formerly  only  total  radiated  energy  was  shown). 

This  still  leaves  continuum  radiation  separate,  however, 
and  in  neither  spectral  bln. 

(B)  scaling  laws  (Eqn's.  1  abc,  2ab)  for  density  and  electron 
temperature-dependence  of  the  low-energy  and  high-energy 
radiative  yields  Y<  and  Y  ; 

(C)  separate  opacity/escape  treatments  for  the  three  spectral 
categories; 


(D)  more  realistic  exchange  of  energy  between  flow  kinetic  energy 
and  thermal,  during  and  just  after  the  "collection"  or  mixing 
process; 

(E)  energy-conserving  averaging  of  Tg  and  when  these  two 
attempt  to  develop  large  rapid  oscillation;  and  finally 

(F)  incorporation  of  a  CRE  radiation  package  for  Argon. 

Yield  studies  for  Aluminum  radiation  vs  wire  array  mass,  initial 
radius,  and  current  penetration  radius  have  been  done,  and  these  comple¬ 
ment  earlier  studies  of  the  variation  with  peak  voltage  and  generator 
risetime.  The  total  yield  is  optimized  at  a  certain  mass,  although  Y> 
optimizes  generally  at  a  considerably  lower  mass,  as  physically  expected 
since  large  masses  cannot  be  heated  to  the  required  temperatures.  Peak 

T  of  course  decreases  with  increasing  load  mass,  while  the  peak  density- 
e 

squared  increases,  until  the  mass  is  too  large  to  be  imploded  in  the 

required  time.  Other  parameters  being  equal,  an  optimum  initial  radius 
TOT  > 

is  seen  for  both  Y  and  Y  ;  larger  initial  radii  have  lower  peak 
temperatures  because  the  assembly  occurs  After  the  current  has  decayed  - 
the  plasma  coasts  in  rather  than  still  being  strongly  squeezed  at  assembly. 
Current-crowbarring  may  weaken  or  shift  this  optimum.  Smaller-radius 
loads  of  the  same  mass  implode  too  quickly  and  tend  to  be  cool,  opaque  and 
subject  to  Simple  Collapse,  with  mostly  blackbody  losses.  Intermediate 
radius  loads  exhibit  well  matched,  Pause-and-Collapse  implosions  in  which 
the  temperatures  are  high  enough  to  forestall  immediate  collapse  and  the 
densities  are  high  enough  for  significant  radiation  but  low  enough  to 
allow  the  radiation  escape  rate  to  match  that  of  kinetic  energy  to  thermal 
conversion. 

The  dynamics  of  the  radiation  pulse  mirrors  the  temperature-dependence 
of  the  radiation  output  curve,  its  multiple  peaks  representing  qualita¬ 
tively  a  time  history  of  Te<  A  single  peak  is  characteristic  of  low 
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temperatures,  ascending  and  descending  the  low-temperature  side  of  the 
L-shell  peak.  Double-peaked  or  triple-peaked  radiation  pulses  tend  to 
be  correlated  with  higher  K-shell  yields,  although  in  real  plasmas 
thermal  gradients,  i.e.,  Tg(r),  may  tend  to  diffuse  and  merge  the  Individual 
peaks.  (SPLAT  is  isothermal  in  the  core  plasma.) 

II.  The  electric  and  magnetic  field  penetration  into  the  Imploding,  time- 
varying  plasma  influences  the  thickness  of  the  current-carrying  region  and 
the  ratio  of  classical  (core)  to  nonclassical  (corona)  currents.  This 
affects  the  circuit  relation,  heating  rates,  and  the  matching  of  the 
implosion  to  the  driver.  To  solve  this  problem,  one  considers  first  the 
coupled  wave-diffusion  equations  in  E  and  B,  and  combines  them  to  form  a 
single  differential  equation  in  a  higher  order  potential  Z,  which  is  a 
generalization  of  the  Hertz  vector  in  classical  electromagnetic  cavity 
theory.  The  equation  for  Z,  its  well-posed  boundary  conditions,  and 
special  techniques  for  its  computational  solution,  have  been  developed 
during  this  past  year’s  contract  effort,  and  are  reported  in  Section  III 
of  this  report,  with  mathematical  details  largely  relegated  to  appendices 
A-D.  As  presently  handled  in  the  SPLAT  code,  magnetic  diffusion  is  frozen 
in  at  an  arbitrarily  chosen  value  at  early  time  (when  the  conductivity 
should  be  low) ,  because  the  plasma  conductivity  rises  to  exclude  further 
field  penetration.  SPLAT  code  results  are  reasonably  sensitive  to  the 
current  penetration,  motivating  us  to  do  this  feature  in  a  more  self- 
consistent  way.  But  the  Hertz  potential  equation  requires  sophisticated 
computer  techniques  which  are  being  developed  under  the  1981  funding. 

III.  Section  IV  of  this  report  treats  beaded  discharges  and  the  implica¬ 
tions  for  the  hard  component  of  the  radiation.  It  assumes  that  the  fully- 
developed  nonlinear  stage  of  the  sausage  instability  necks  the  plasma 
column  off  in  many  places  to  produce  relatively  low-density  regions. 

The  condition  for  nonclassical  conduction  in  these  regions  is  remarkably 
lenient  at  high  temperature,  so  we  Investigate  the  consequences  of  the 
regions  being  nonclassical  at  all  radii.  At  high  E-flelds,  this  gives 
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lower  current  than  with  classical  conduction,  but  gives  greatly  enhanced 
ohmic  heating,  to  tens  of  keV,  and  this  energy  is  deposited  in  the 
denser  beads,  producing  hard  radiation.  Depending  on  the  E-field  pene¬ 
tration,  which  is  enhanced  by  the  unstable  fluctuations,  this  effect  may 
produce  high  power  levels  of  such  radiation,  and  offers  the  prospect  of 
copious  K-shell  radiation  from  plasmas  whose  mass  is  too  large  to  allow 
the  temperatures  required  for  thermal  K-shell  emission.  The  efficiency 
is  high  if  much  of  the  voltage  drop  occurs  in  the  low-density  regions 
between  beads,  and  if  these  regions  are  long  enough  to  allow  adequate 
heating  before  deposition. 

A  short  description  of  ongoing  work  is  contained  in  Section  V  of 
this  report.  Two  areas  are  under  development:  high-order  interpolation 
techniques,  which  are  necessary  for  computer  solution  of  the  magnetic 
field  diffusion  problem  but  also  useful  in  radiation  modeling  and  MHD 
codes,  and  the  modification  of  a  ID  MHD  code  to  treat  the  dynamics  of 
low-mass  higher-temperature  imploding  plasmas  which  are  not  "strongly- 
coupled"  and  thus  not  properly  described  by  the  SPLAT  code.  (Typical 
Argon  implosions,  for  example,  probably  do  not  have  the  viscosity  neces¬ 
sary  for  the  SPLAT  description.)  Because  of  the  potential  importance 
of  the  field  fiffusion  and  radiation  processes,  we  are  working  toward 
making  the  MHD  code  interface  easily  with  the  generalized  Hertz-potential 
solver  and  the  local  temperature  equations  with  CRE  radiation. 
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CHAPTER  II.  SYSTEMATIC  YIELD  TRENDS  RESULTING  FROM 
THE  PRESENT  CORE-CORONA  MODEL 

A.  Code  improvements  over  this  contract  period 

The  existing  simplified  core— corona  implosion  model2,  has  been  modified 

to  provide  spectral  resolution  in  the  radiative  yields  computed  during  the 

course  of  a  single  implosion.  Extending  previous  techniques,  the  radiative 

emission  from  an  aluminum  plasma  of  fixed  density  and  radius  is  calculated 

as  a  function  of  electron  temperature  using  the  collisional-radiatlve  equili- 

2 

brium  (CRE)  model  due  to  Dus ton  and  Davis  .  The  size  of  the  cylindrical 
radiator  is  set  to  a  small  value,  e.g.  0.05  cm,  and  the  output  radiated  energy 
decomposed  into  three  categories:  "L-shell"  or  hv  <  1  keV,  "E-shell”  or 
hv  >  1  keV,  and  recombination  continuum.  One  may  then  derive  an  equivalent 


volume  emission  rate  (e)  from  this 

represented  by 

single  line  source  and  these  results  are 

\ 

e<(n°,  T  )  -  exp  j  S 

L  e  Ik 

vl 

in(Te/0.3)]k| 

(la) 

e>(n°,  Tg)  -  exp  j 

V 

[  ta(Te  /0.7)]k| 

(lb) 

eFB(n°,  Tfi)  -  exp  j]C 

.  FBI 
bk 

[in(Te/  0.5)] k| 

(lc) 

at  the  reference  density  n^.  This  emission  is  then  subjected  to  a  density 
scaling  appropriate  to  each  radiation  category,  in  particular,  a  reasonable 

3 

approximation  is  obtained  using 

€>,FB  ^  (nj/1018)1’5 

nI 

>  1018 

(2a) 

(nj/1018)2 

"I 

<1018 

e<  ~  (nj/1016)1,5 

nl 

>  io16 

(2b) 

^  (nj/1016)2 

“I 

<  1016  . 
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Plots  of  the  spectral  representations  e 


>  FB 

e  ,  £  are  shown  in  Figure  1 


as  P\  P>, 


.FB 


P  respectively  arising  from  a  CRE  calculation  with 


19  -3 

1.0*10  cm  and  a  radius  of  0.05  cm. 

FB 


The  curve  P  ,  is  the  sum  of 
rad  - 

the  Individual  representations  P' ,  and  can  be  seen  to  reproduce  the 

original  CRE  results  (marked  by  a  "q")  to  within  5  Z  or  better  over  the 

temperature  domain  of  interest  here.  The  equivalent  volume  emission  rate 

A  2  , 

from  a  black  body  of  the  same  size  (a  Tg  /  w  r q  ;  is  plotted  for  comparison. 

Using  this  c(n^,  T  )  as  a  source  term,  the  net  radiative  loss  from  any 

plasma  density  profile  |n^.(r^,t)  and  T^(t)J  generated  bv  the  core-corona 

model  can  be  developed  once  the  proper  opacity  corrections  are  calculated 

in  a  local  approximation.  At  selected  locations  r^  within  the  density 

profile  e<(nI(ri)  T^) ,  e>(nI(ri>  T^) ,  eFB(nI(ri> ,  T^)  and  e^OijCr^  Tx) 

(Bremsstrahlung)  are  calculated  and  subjected  to  an  additional  attenuation 

based  on  the  probability  of  escape  (P  )  from  the  specified  location  r. . 

e  14 

The  P£  depends  upon  the  optical  depth  T  ^  and  an  accurate  approximation 
to  this  escape  probability  is  to  compute  optical  depths  along  each  of  two 
paths,  inclined  to  the  radius  vector  along  a  chord  representing  the  mean 
angle  of  emission  for  a  cylindrically  symmetric  radiator,  and  average  the 
P^  over  the  two  path  lengths.  This  implies  that  an  effective  optical  depth 
can  be  assigned  to  r^  once  the  line  Integral 


r* */: 


d& 


TVjd) 

ni(Ti> 


is  computed  over  the  short (S)  and  long(L)  path  to  the  plasma  surface.  Each 
point  r^  is  thus  assigned  two  optical  depths  for  each  spectral  category, 
e.g.  t*  _  and  t*  ,  with  these  optical  depths  computed  from  t.  l.  T. 

i>2>  X|L  1»*> 

The  appropriate  scattering  cross-sections  for  a  0.5  keV  "L-shell"  line, 

a  1.5  keV  "K-shell"  line,  Al  photo-ionization/recombination,  and 

<  >  FB  FF 

Bremsstrahlung  provide  Tg  ^r^,  Tg  L(r1>,x  /g  L(ri>  and  *s  ^(rj)  respectively 
once  l.  c  and  l.  .  are  known. 

The  final  radiative  loss  from  each  point  in  the  density  profile  can  be 
written 
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Pf#d  (erg/cm3  sec)  as  a  function  of  T#  (keV) 

{from  a  0.05  cm  radius  cylindrical  Al  plasma  with  n.  ■  1.0-1019  cm*3) 


I  I  I  1  1  1  iJ _ I _ I _ I  I  I  I  I 

13  .04  06  .08  .10  .2  .3  .4  .6  .8 
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Figure  1 


FF  FF 


-  P  •  Y(e’TS  +  e'^  )  +  P 


(3) 


+  e>  •  2(Pe(TS>)  +  Pe<Ti>  +  C<  *  7(Pe(TS<>  +  Pe(TL<)) 


and  Chls  function  is  radially  integrated  to  produce  the  radiant  energy  loss 

at  each  time  step.  As  a  spectral  diagnostic,  the  radially  integrated  power 
^  v  FF  FB 

loss  for  lines  (P  ,  F  )  and  continuum  (P  +  P  )  is  reported  separately 
and  time  integrated  as  well  to  provide  accumulated  yields  in  the  three 
categories  "L-shell"  (h^  <  1  keV) ,  "K-shell"  (hy  >  1  keV)  and  all  continuum 
(FF  +  FB) . 

Within  this  local  approximation  the  model  plasma  is  represented 
(radiatively)  as  a  superposition  of  uncoupled  CRE  line  sources  attenuated 
by  the  local  probability  of  escape  appropriate  to  each  spectral  category. 

This  P  (r^)  is  calculated  self-consistently  with  respect  to  the  density  and 
temperature  profile  generated  by  the  MHD  plasma  evolution,  and  the  individual 
line  source  strength  e(nI>  T  )  is  assigned  self-consistently  with  respect  to 
this  evolution.  In  this  context  the  local  radiation  approximation  is  in 
keeping  with  the  strongly  coupled  fluid  limit  used  within  the  core-corona 
model  (at  the  outset  of  this  research)  to  simplify  the  MHD  plasma  evolution. 

In  a  local-radiation  limit  one  has  the  assumption  that  the  radiative  diffusion 
of  energy  is  so  rapid  that  all  the  superposed  CRE  sources  are  in  equilibrium 
with  each  other  so  that  there  is  no  need  to  couple  them  and  Iterate  in  order 
to  establish  the  excited  state  populations  (note  that  nj(r,t)  must  be  a 
conserved  quantity  in  that  iteration).  In  fact,  for  optically  thick,  highly 
collislonal  systems  it  is  this  very  rapid  radiative  diffusion  which  helps 
support  the  isothermal  plasma  temperature  profile.  In  short,  the  thicker 
the  plasma  the  more  appropriate  this  local  radiation  limit  becomes,  and  the 
method  can  be  extended  to  an  arbitrary  nwber  of  spectral  categories. 

A  similar  development  has  been  carried  through  for  Ar  and  the  SPLAT 
code  is  now  equipped  with  a  radiation  package  which  models  either  A1  or 
Ar  depending  on  the  input  choices.  The  equivalent  CRE  source  profile 
for  Ar  is  shown  in  Figure  2. 
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Pra<j(erg/cm3  sec)  as  a  Function  of  Te(keV) 
from  a  0.001  cm  Ar  plasma  with  ■  1.0‘IQ1*  cm 
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A  second  area  of  improvement  is  in  the  thermalization  of  the  kinetic 
energy  in  the  flow  field  V*)  during  the  collection  process.^  The 

present  technique  involves  an  explicit  calculation  of  the  energy  residing 
in  this  flow  at  each  time  step  during  collection  and,  when  the  flow 
energy  at  the  new  time  step  is  less  than  that  at  the  old,  the  energy  dif¬ 
ference  is  thermalized  either  by  setting  a  heating  rate  AE/At  or  by 
instantaneously  boosting  the  core  internal  energy  by  AE.  In  addition, 
energy  deficits  at  the  end  of  the  collection  process  are  thermalized  over 
a  time  scale  At  -  (a-r  )/(c,-a)  characteristic  of  a  shock  transit.  This 
models  the  final  momentum  propagation  to  the  outer  regions  of  the  plasma, 
as  required  to  bring  the  velocity  profile  into  a  centrally  stagnant  form 
when  the  deceleration  of  the  annulus  proceeds  from  the  center  outward. 

A  final  improvement  is  the  installation  of  an  internal  energy 
conserving  average  for  Tg  and  T^  in  the  core  when  the  energy  exchange 
time  is  very  short  and  one  seeks  to  stabilize  temperature  oscilla¬ 

tions.  The  quantity  u  *  T^  +  Z(Tg)Te  is  now  kept  fixed  whenever  rapid 
exchange  oscillations  imply  the  need  to  average  Tg  and  T^.  The  (usually 
slow  but  sometimes  significant)  electron-ion  heat  exchanges  are  now  also 
included  in  the  coronal  temperature  calculation  (zT2) . 

B.  The  Standard  Yield  Experiment 

In  order  to  mimic  a  real  experiment  with  our  present  core-corona 
model  several  constraints  must  be  kept  in  mind.  First,  there  is  every 
reason  to  believe  that  present  experiments  are  most  surely  not  cylindrical 
in  their  implosion.  The  wire/puff  gas  plasmas  tend  to  close  on  the  axis 
like  a  zipper  due  to  inhomogeneities  in  the  electromagnetic  fields  within 
the  diode  cavity.  Since  no  axial  variations  are  followed  in  the  present 
code,  one  finds  that  accurate  models  of  machine  driving  pulses  provide  an 
excess  of  energy  to  the  1-D  system  and  assemble  it  violently  on  axis. 

Under  some  conditions  this  results  in  many  bounces  of  the  model  plasma, 
which  are  unrealistic, and  the  yield  figure  extracted  from  the  model  must 
take  this  into  account.  Second,  the  present  model  does  not  treat  the 
late  time  behavior  of  the  discharge  when  "sausage"  and  "kink"  instabilities 
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are  often  apparently  present.  Since  the  sausage  mode,  in  particular, 
probably  has  a  significant  impact  on  the  (time  dependent)  discharge 
Impedance  and  thus  on  its  radial  confinement  and  compression,  the  yield 
figure  extracted  from  the  present  code  must  avoid  the  post-implosion 
phase.  Third,  the  present  model  does  not  attempt  to  treat  the  "crowbar" 
effect  active  in  many  existing  diodes,  which  keeps  the  current  peaked 
at  a  large  value  while  the  trapped  magnetic  energy  decays.  In  a  device 
that  does  not  short  out,  one  expects  the  current  to  decay  smoothly  as  the 
machine's  driving  voltage  decays  and  this  behavior  is  indeed  obtained 
from  our  circuit  equation.  The  weakening  of  the  current  allows  the 
present  model  plasma  to  expand  more  easily  than  one  might  expect  in  the 
real  experiment  and  again  renders  the  yield  figure  developed  over  the 
later  phases  of  the  implosion  inapplicable  to  existing  experiments. 

A  simple  yield  observation  can  be  defined  that  avoids  these  caveats 
if  one  simply  presumes  that  some  time  window  dependent  on  the  bulk  plasma 
motion  can  be  extracted  experimentally.  For  purposes  of  a  general  comparison 
let  Yg  be  defined  as  that  yield  (in  any  spectral  category)  developed  through 
a  single  compression  and  subsequent  expansion  of  the  load,  i.e. 

/e 

V,  4|  dt  P  (t.)  (5) 

E  J  i  rad 

o 

with  t£  the  (expansion)  time  when  the  radially  outward  flow  velocity  first 
falls  to  zero  after  the  Initial  compression.  This  choice  of  a  yield  measure¬ 
ment  removes  all  late  time  behavior  and  represents  a  realizable  measurement 
once  the  experimental  radiative  output  is  correlated  in  time  with  some 
measure  of  the  plasma  flow  velocity.  The  present  model  can  thus  be  used  to 
study  Y£  ,  Yg  ,  Yg  as  functions  of  some  20  independent  model  input 
parameters. 

Our  recent  work  focuses  on  only  two  of  these  input  parameters,  but  ones 
which  are  clearly  related  to  easily  adjusted  experimental  initial  conditions. 
These  are  the  total  mass  of  the  load  (Mq)  and  the  location  of  the  density  maxi¬ 
mum  (r°),  upon  transition  to  a  plasma  state,  and  are  easily  varied  as  parameters 

IB 

of  the  initially  self-similar  plasma  model.  These  are  experimentally  equiva¬ 
lent  to  the  number  of  wires  (of  fixed  diameter)  in  a  load  array  and  the  initial 
radial  location  of  the  load  array.  In  order  to  complete  the  specification  of 
these  studies,  a  list  of  the  other  relevant  input  parameter  ranges  is  shown 
in  Table  I. 


Table  I.  Fixed  Ittitial  Parameters 


Diode  Inductance, 

13.5  nh 

Generat  or  Impedance , 

0.5  fl 

Maximum  V  (t) 
oc 

♦  4.5  MV 

Time  to  reach  V  Bax 
oc 

150  ns 

Initial  current  penetration  parameter*  g  as 20  — 

—  5 

g  —  (r^a)2  /  l-trj/a)2 

Initial  plasma  temperature 

15  eV 

Initial  annulus  width 

*  0.14  cm 

Return  current  radius 

1 _ 

6  cm 

*this  adjusts  the  tine-average  core/corona  current  partition 
0(t),  with  larger  g  values  producing  smaller  0  values. 


In  modeling  experiments  which  vary  the  r°  and  M_  parameters,  the  total 

m  o 

number  of  ions  should  correspond  to  the  mass  of  a  typical  load,  and  masses 

corresponding  to  [6  -*■  24]  1.5  mil  A1  wires  3  cm  in  length  were  used.  This 
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provides  a  (1.125  -*■  4.50)  x  10  range  in  the  total  ion  number  for  the 
load  or,  equivalent ly,  a  load  mass  range  of  [168  -*■  673  yg  J  cm].  The  initial 
profiles  containing  a  fixed  number  of  ions  are  chosed  so  as  to  (i)  con¬ 
fine  the  ions  within  a  0.2  cm  (or  less)  thick  annulus,  (ii)  provide  for 

surface  velocities  on  the  order  of  the  sound  speed  at  T,  (t  )  ~ 15  eV, 

i,e  o 

and  (ill)  assign  surface  accelerations  which  are  compatible  with  the 
pressure  gradients  and  initial  Jg  x  Bq  stress  profiles.  Second,  the  peak 
density  of  the  profile  is  chosen  to  correspond  with  the  Initial  load 
radius  and  has  been  varied  over  [0.3  -*■  2.0]  (cm),  covering  popular 
experimental  choices. 


A  third  Important  parameter  in  the  profile  is  rj(<  a) »  the  effective 
skin  depth  of  the  plasma.  This  choice  cannot  be  set  easily  by  experiment 

t 

but  it  is  a  very  Important  parameter  in  the  yield  studies.  A  range  in 


_  of  T5.0  -►  20]  has  been  examined  in  recent  work.  The  wide 

g  - -  b 

l-(rl®)2  \ 


domain  for  g  is  necessary  because  so  little  is  known  as  to  its  proper  value. 
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In  summary  che  available  radius  and  vire/mass  combinations  are 
denoted  in  later  figures  by 


rm  “X 

wires  (M  ) 
o 

6 

8 

12 

16 

24 

2.05 

A6 

A8 

A12 

A16 

A24 

1.58 

B6 

B8 

B24 

0.97 

C6 

C8 

C12 

C16 

C24 

0.68 

D6 

D8 

D12 

D16 

D24 

0.30 

E6 

E8 

E12 

E16 

E24 

with  Vp  and  g  continuously  variable. 

C.  Dynamics  of  Typical  Implosions 

As  a  general  rule  the  core/corona  implosion  falls  into  one  of 

three  distinct  trajectory  types.  These  are  illustrated  below  in  Figures 

3,  4,  S  for  the  case  codes  and  legends  shown.  The  demarcations  between 

the  SC/PC/BC  trajectories  in  the  experimental  domain  [r°,  M  ,  V  ,  g]  are 

m  o  p 

rather  complex  as  the  radius  decreases  and  not  yet  mapped  out.  For  the 
larger  radii,  however,  the  transition  through  PC  to  rather  clear  BC 
trajectories  is  progressive  as  the  load  mass  is  decreased. 

The  general  importance  of  the  corona  is  also  a  strong  function  of 

radius  r°  and  g.  The  smaller  values  of  r°  tend  to  produce  significant 
m  m 

current  transport  in  the  corona  at  peak  compression  for  moderate  g  values. 

The  larger  values  of  r°  tend  to  allow  most  of  the  current  to  flow  in  the 

n 

core  most  of  the  time  unless  the  g  values  are  very  large  ( 40  -  80) . 

Clearly  the  details  of  this  dependence  are  of  interest  in  load  design  but 

more  initial  condition  points  in  r°  must  be  developed  for  such  a  study. 

m 

The  physical  origin  of  the  effect  is  rooted  primarily  in  two  things. 

First, small  r°  loads  tend  to  compress  more  and  thus  drive  g  up  as  the 

heating  due  to  assembly  occurs.  Second,  the  hot  dense  core  tends  to 

release  energy  into  the  corona  when  its  surface  density  gradients  are  so 

stiffened.  These  effects  in  concert  produce  a  drop  in  6  =  I  _  JL  . 

■'  cor*/  coux 
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Figure 


tory  (Case  C12,  g  *  13.1) 


At  large  r  values  the  assembly  heating  generally  prevents  stiff  gradients 
m 

by  holding  off  the  compression.  These  gradients  are  available  later  when 
the  core  plasma  has  been  cooled  radiatively  but  then  the  core  has  little 
energy  to  release  Into  the  corona.  Hence  these  effects  when  separated 
produce  little  decay  In  8. 

The  radiation  pulse-width  In  these  Implosions  can  be  seen  to  vary 
over  a  fairly  wide  range  5-30  nsec.  The  shorter  pulses  generally  are 
obtained  from  cooler,  continuum  emitting  plasmas  while  some  of  the  larger 
pulses  obtain  when  the  timing  of  assembly  heating,  compression  and  peak 
current  all  conspire  to  produce  a  hot  plasma  [600  eV  -  1500  ev]  that  is 
momentarily  confined  or  quickly  recompressed  (in  a  PC  sense)  by  large 
currents  confined  at  small  radius.  For  the  standard  voltage  rise  time 
mentioned  above,  only  the  lower  values  of  r°  (<  1.0)  can  give  rise  to  this 
sort  of  pulse  because  the  load  must  be  able  to  drift  into  a  nearly  assembled 
state  just  before  the  peak  voltage  and  current  rise  to  heat  and  trap  it. 

D.  Standard  Yield  Experiment  Parameter  Studies 

Perhaps  the  easiest  variation  to  examine  is  the  dependence  of 
yield,  Yc  in  various  spectral  categories,  on  the  generator  peak  voltage. 

In  Figure  6,  this  quantity  is  plotted  for  values  of  <  3.0  MV  using 
initial  case  C12  with  an  assumed  g  “  9.0  for  each  of  three  spectral 
categories:  Yg°tal,  Yg,  Yg.  Also  shown  are  the  peak  temperatures 
achieved,  T^  and  the  peak  compression,  n^,  over  the  time  frame  of  the 
Implosion.  The  generator  voltage  must  be  sufficiently  large  to  "ignite" 
the  load,  but  once  this  threshold  is  achieved  further  voltage  increases 
do  not  improve  the  yield  significantly.  This  results  from  a  very  much 
weaker  compression  for  higher  voltages  (which  give  higher  implosion  speeds) 
as  evidenced  by  a  smooth  increase  in  peak  temperature  and  drop  in  peak 
density. 

A  second  fundamental  variation  is  that  of  dependence  on  the  array  radius 
r°.  In  Figure  7,  the  sequence  of  initial  cases  E12,  D12,  C12,  B12,  A12 
has  been  used  to  produce  the  same  information  as  above:  YgOCa^,  Yg,  Yg, 
T1,ni,  at  a  fixed  voltage  VQ  ■  1.5  MV  and  g  ■  9.0.  Here  there  is  a  clear 


.3  -- 


optimum  in  r°  because  of  a  combination  of  kinematic  and  dynamic  effects, 
m 

At  large  r°,  even  though  S  rises  rapidly  to  near  0.97  from  a  small  initial 
m 

value,  the  load  has  too  large  a  distance  to  cover  in  the  fixed  rise  and 

fall  time  of  the  driving  voltage.  This  produces  the  larger  currents  while 

the  conductor  is  at  larger  radius  and  hence  does  not  give  the  larger  values 

of  J  x  Bn  that  can  be  obtained  with  smaller  r°.  Bringing  r°  down  thus 
z  a  m  m 

gives  a  "kinematic"  optimization  and  improves  the  yield,  peak  temperature 
and  compression.  On  the  other  hand  if  r°  is  too  small  the  load  experiences 
strong  inward  forces  but  does  not  have  enough  distance  of  travel  to  allow 
these  forces  to  deliver  much  kinetic  energy.  This  is  a  dynamic  degrada¬ 
tion  of  the  implosion  quality,  caused  by  being  too  close  to  the  bottom  of 
the  effective  potential  well  created  by  the  generator. 

Finally  the  variation  with  load  mass  is  shown  in  Figure  8,  generated 
by  initial  cases  C6,  C8,  C12,  C16,  C24  at  fixed  voltage  Vq  and  g  =  9.0.  Also  for 
comparison  this  mass  variation  is  shown  in  Figure  9  with  r^  fixed  rather  than 
g.  In  this  later  sequence  g  will  change  slightly  at  tQ  because  the  profiles 
which  evelopea  smaller  mass  do  not  spread  as  far  from  the  peak  r  °  as  those 

ID 

of  larger  mass.  Physically  the  sequence  with  fixed  g  correpsonds  to  early 
histories  which  imbed  the  current  in  a  fixed  fraction  of  the  load  volume,  vnlle 
the  sequence  with  fixed  r j  corresponds  to  early  historic-,  that  preface  thinner 
skin  depths  on  smaller  masses  due  to  a  more  rapid  heating  <?hich  freezes  the 
current  diffusion  before  it  has  a  chance  to  penetrate.  Dynamically  the  yield 
will  degrade  for  both  large  and  small  masses.  Large  masses  cannot  be  accel- 
rated  well  and  thus  build  no  kinetic  energy  reservoir  with  which  to  drive  the 
radiation.  Small  masses  heat  too  rapidly  and  hold  off  the  compression  of  the 
load,  degrading  the  yield  because  too  low  a  density  is  achieved. 

Finally,  in  Figure  10,  the  variation  of  the  yield  with  g  is  shown  in  order 
to  illustrate  the  general  sensitivty  of  the  implosion  quality  to  the  detailed 
microscopic  assumptions  underlying  the  calculation  of  the  load  dynamics. 

Larger  g  values  produce  more  rapid  compressions  and  higher  peak  temperatures 
because  they  create  higher  surface  stresses;  roughly  similar  total  currents 
are  confined  to  a  narrower  annulus,  increasing  J  x  B. .  For  very  large  g 

Z  o 

values,  of  course,  the  core  current  will  decay  (8  -*•  0)  because  infinite 
current  densities  are  not  allowed  due  to  coronal  growth. 
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Yield  vs.  Mass,  fixed 


Figure  10.  Yield  vs.  g(t 


E.  Testing  the  Yield  Dependence  on  a  Possible  Scaling  Parameter 

Recent  Maxwell  Laboratory  experiments^  have  examined  the  utility 
_  o  ^ 

of  the  parameter  =  r^  Mq  as  a  measure  of  implosion  performance.  It  has 
been  shown  to  be  a  poor  indicator  of  implosion  quality  with  yields  decay¬ 
ing  rapidly  and  peak  temperatures  increasing  as  one  progresses  to  larger 
values  of  r°  at  constant  K^.  The  load  masses  in  the  Maxwell  experiments 
are  smaller  than  those  used  in  the  studies  described  above,  but  SPLAT 
calculations  (at  the  larger  masses  used  before)  indicate  that  fails 
as  a  measure  of  implosion  performance  for  the  model  plasma  also. 

Using  the  standard  generator  impedance  and  inductance  from  Table  1 

and  adjusting  the  driving  voltage  waveform  model  for  a  80  ■+■  150  nsec 

rise  time  to  a  peak  voltage  of  1.5  -*■  3.0  MV,  a  sequence  of  SPLAT  calcula- 
o 

tions  at  variable  with  fixed  —  17.8  £cm.yg]  shows  a  trend  similar 

to  the  experimental  behavior.  In  Figure  11  these  results  are  shown  as  a 

function  of  array  radius.  The  yield  falls  because  of  a  systematic  trend 

toward  overheating  at  large  radii,  if  sufficient  mass  is  removed  to  hold 

K.^  constant.  The  greater  temperatures  reverse  the  implosion  too  soon  and 

prevent  compression  to  densities  sufficient  for  appreciable  radiative  yield. 

According  to  the  SPLAT  code  dynamics,  the  proper  contour  of  constant  yield 

would  involve  a  slower  power  of  mass,  K  —  r°  m\  with  y  <  \  in  order  to 

m  o 

prevent  the  overheating. 

This  failure  of  K  to  represent  a  figure  of  merit  corresponding  to 

T  >  1  7 

yield  Y  or  Y  is  due  probably  to  the  origin  in  an  optimization  based  on 

constant  levels  of  power  coupling  from  the  generator  to  the  load.  Since 

the  input  power  levels  are  only  weakly  related  to  the  time  integrated 

radiative  loss,  or  yield,  this  parameter  is  suspect  at  the  outset. 
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CHAPTER  III.  ELECTRODYNAMIC  MODELS  OF  ANNULAR  PLASMA  LOADED 
DIODES  BASED  ON  A  GENERALIZED  HERTZ  VECTOR  POTENTIAL 


A.  Introduction 

The  common  use  of  a  circuit  equation  (usually  derived  from  Faraday's 
law)  in  MHD  models  of  imploding  annular  plasma  loads  is  grounded  in  very 
serious  simplifications  of  the  actual  electrodynamics  in  pulse  driven 
diodes.  The  most  important  simplifications  involve  the  use  of  an  essen¬ 
tially  electrostatic  approximation  for  the  electric  field,  the  use  of 
magnetic  field  diffusion  theory  in  the  presence  of  an  inhomogeneous,  time 
varying,  convecting  conductor  (the  plasma  load),  and  the  neglect  of  the 
tensor  character  of  the  collisional  plasma  conductivity  due  to  the  imbedded 
azimuthal  magnetic  field.  Moreover,  the  geometry  of  actual  machines  is 
complex  enough  to  preclude  detailed  electromagnetic  theoretical  analysis 
in  anything  less  than  three  spatial  dimensions. 

Many  difficulties  can  be  removed  from  the  conventional  theoretical 

g 

treatments  by  simply  insisting  on  a  more  rational  diode  geometry.  A 
simpler,  more  z-symmetric  diode  feed  would  obviously  simplify  analysis 
and  prediction  of  the  experimental  results.  Significantly,  the  choice  of 
new  diode  geometries  may  be  guided  by  theoretical  analysis  using  a  generali- 
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zation  of  the  Hertz  potentials  ’  which  are  intrinsically  suited  to  the 
problem  of  fluid  conducting  plasma  in  a  wave-guide.  The  choice  of  boundary 
conditions  can  select  particularly  simple  electromagnetic  field  solutions 
and  load  configurations  that  admit  a  relatively  tractable  analysis  over 
a  rather  general  class  of  waveguide  geometries  and  electrodynamic  load 
properties. 

The  new  field  equations  for  the  generalized  Hertz  potentials  are 
equivalent  to  the  Maxwell  equations  which  give  rise  to  them,  but  they  are 
simpler  in  that  a  single  vector  field  is  derived  that  represents  E  and  B. 

The  solution  of  a  single  wave  equation  (or  wave-diffusion  equation)  is 
then  sufficient  to  provide  all  the  electromagnetic  field  components.  This 
new  vector  field  Z(x,t)  is  equivalent  to  the  usual  Hertz  vector  potential 
only  in  source  free  regions  and  therefore  represents  a  non-trivial  generali¬ 
zation  of  the  original  field.  The  choice  of  Z(x,t)  is  in  fact  the 
generalization  of  the  familiar  process  of  making  a  convenient  gauge  choice, 
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allowed  by  the  arbitrary  lamellar  component  of  the  ordinary  vector  poten¬ 
tial  A.  Here  Z_  will  be  chosen  so  as  to  remove  the  need  for  two  inhomo¬ 
geneous  Maxwell  equations,  and  to  replace  them  with  a  single  equation  to 
determine  Z. 

B.  The  fundamentals  Revisited 

Such  a  program  of  reduction  can  always  be  carried  out  when  the 
field  sources  are  prescribed  functions,  i.e.,  p  =  p(x,t) ,  .J  =  >J(x,t)  are 
given  and  the  fields  external  to  these  sources  are  required.  In  this 
case  one  must  seek  a  vector  potential  £(x,t)  which 

(i)  allows  the  definition  0(x,t)  «  -  •  Z  (x,t)  to  represent 

a  gauge  choice  relating  0,  A  (x,  t)  (the  usual  vector 
potential)  and, 

(ii)  reduces  the  number  of  inhomogeneous  field  equations  from 
two  (for  0  and  A  separately)  to  one  for  the  new 
potential  Z. 

Performing  the  usual  elimination  of  the  homogeneous  Maxwell  equations  by 
introducing  0  and  A,  the  two  remaining  relations  become 

£.(\^$  +  c^3tA)"-4ir  p  (4a) 

V  x  V  x  A  +  c’1  3t  (V  4  +  c"1  3t  A)  -  +  4  ir  c-1  J.  (4b) 

Now  choose  the  generalized  Hertz  vector  Z  so  that 
d  -1 

A  -  c  1  3t  Z  (5) 

d  -1 
and  if  0  =  -  V^Z ,  then  a  familiar  gauge  choice  (V*A  +  c  3  0  *  0)  is  implied. 

If  one  now  replaces  0  and  A  in  (4a,  b)  with  the  equivalent  Z  terms,  the 

equations  become 
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(6a) 


□  (7-Z)  -  4  ir  p 
3  □  Z  *  —  4  it  J. 

1  ~  (6b) 

The  procedure  for  choosing  Z  is  formalized  in  Appendix  A,  which  points 
out  the  role  of  a  Hertz  vector  potential  in  the  conventional  electro¬ 
dynamics  . 

Let  the  source  field  J  be  represented  by  the  time  derivative  of  the 
charge  transferred  per  unit  area  in  any  spatial  direction  over  the  entire 
time  domain  of  interest,  viz. 


J.  ■  T  (x,  t) ,  or 
T  |  f  dt  J(x,  t)  +  T°  (x)  . 


Using  the  continuity  relation,  it  is  easily  shown  that  (x,  t)  reduces  the 
problem  to  a  single  wave  equation 

Q2  Z  -  -  4  ir  T,  (8) 

because  (6a)  is  subsumed  in  the  divergence  of  (8)  so  long  as  T° (x)  is  a 
solenoidal  field.  In  terms  of  the  potential  Z  and  transfer  vector  T,  the 
observable  fields  are  written 

E  m  2.  x  E.  x  —  (*»  -  4ir  ^ 

B  -  c"1  7  x  3t  Z  (x,  t)  (9b) 

and  the  uniqueness  of  £  (x,  t)  is  guaranteed  because  the  specification  of 
p  or  i.e.  ]T  (x,  t) ,  and  the  initial  fields  £  (x,  tQ) ,  £  (x,  tQ)  at  some  time 

tQ,  can  be  seen  to  imply  Cauchy  open-surface  boundary  conditions  for  the  wave 
equation  (8).  The  uniqueness  of  Z  is  also  discussed  in  Appendix  A. 


In  a  waveguide  formed  by  two  Infinite  plane  parallel  perfectly  conducting 
surfaces,  the  Hertz  vector  Z  represents  a  variety  of  electromagnetic  wave 
solutions.  However,  given  the  azimuthal  symmetry  and  perfectly  conducting 
guide  surfaces,  only  a  particular  subset  are  of  interest  in  analyzing  and 
describing  the  evolution  of  an  annular  plasma  load.  One  seeks  radially 
varying  TEM  modes  which  limit  asymptotically  as  r  -♦■<*>  to  a  superposition  of 
incoming  and  outgoing  cylindrical  waves  with  (perhaps)  a  quasi-static  component 
added.  At  large  radius,  the  vacuum/waveguide  boundary  conditions  therefore 
restrict  consideration  to  those  solutions 

Z  -  z  J  cC  (r,  t)  +  5*  (z,  t) j  (10) 

where  t,  A  denote  the  solenoidal  and  lamellar  portions  of  the  Hertz  vector. 

The  radial  component  of  Z_  must  vanish  because  if  one  insists  that  the  scalar 
potential  (-  V_  •  Z)  exhibit  no  radial  gradients  (in  consonance  with  the  boundary 
conditions) ,  then  the  E  arising  from  the  solenoidal  portion  of  Z  can  vanish  if 
and  only  if  Z^  itself  is  everywhere  zero. 

The  delicate  question  is  how  to  preserve  this  solution  as  it  collides  with 
and  is  reflected  by  a  convectlng  quasineutral  (p  m  0)  plasma  load  which  exhibits 
strong  conductivity  gradients.  The  plasma  must  possess  and  self-conaistently 
retain  several  symmetries  in  order  to  insure  retention  of  Z  in  the  form  given 
above.  The  central  requirement  is  that  E,  be  orthogonal  to  Vo.  Since  JE  arises 
from  the  independent  fields  0  and  A,  it  is  sufficient  to  require  that  0  have 
no  radial  gradients  and  that 


V  X  V  -  0, 

(11a) 

V  x  v_  O  -  0  , 

(lib) 

o  V  •  V  x  V  x  A  -  0, 

(11c) 

where  one  has  in  the  rest  frame  of  the  convectlng  plasma  (superscript  ') 

}• 

orthogonal  to  the  gradients  of  cr(x,t)  =  a(x,t)  y(V(x,t))  and  to  the 
convection  field  V(x,t). 


4 


lx 


o  (x,  t)  Y  (V  (x,  t)) 


|  -  V  0  - 


-1 


V 

3_  A  +  — 
t  —  c 


V  x  A 
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The  next  step  is  to  try  to  carry  out  the  program  used  above  when  p 
and  J  were  given  functions.  In  this  case  the  choice  of  Z(x,t)  is  not  so 
transparent;  one  must  construct  the  proper  functional  A(Z)  which  provides 
the  results  (i)  =  -  V*Z^  and  (ii)  only  one  (in  this  case  homogeneous) 

propagation  equation  is  non-trivial.  Such  a  construction  can  be  done  using 
the  Green's  function  solution  of 

V*  G(x,x)  =  -  4tt<5(x-x) 


which  provides  the  proper  boundary  conditions  in  the  waveguide.  For  the 
plane  parallel  perfect  conductors  and  the  axially  independent  cylindrical 
wave  solution  at  large  radius,  this  requires  G(x,x)  =  G(r,r,  0-&) ,  a 
single  line  source  located  at  (r,$);  but  the  method  outlined  in  Appendix  B 
is  applicable  in  any  geometry. 

The  solution  algorithm  for  A(Z )  and  Z_  (as  discussed  in  Appendix  B) , 
uses  the  mathematics  of  space-ordered  exponential  operators,  and  appears  to 
provide  a  relatively  straightforward  computational  means  of  solving  the 
otherwise-intractable  problem  of  fields  in  the  presence  of  a  convecting 
inhomogeneous  plasma  load  with  time-varying  conductivity.  As  long  as  the 
diode  surfaces  are  approximated  as  perfectly  conducting,  the  solution  demands 

(TEM)  modes  governed  by  a  single  transverse  component  of  the  vector  Z,  i.e. 

I 

C  -*■  0  everywhere.  When  the  perfectly  conducting  waveguide  and  the  plasma 
are  configured  such  that  (a)  radial  plasma  convection  currents  are  divergence¬ 
less  (i.e.,  plasma  waves  in  the  ambipolar  space  charge  can  be  neglected), 

(b)  the  plasma  gradients  are  radial,  and  (c)  the  plasma  is  axially  uniform, 
then  all  axial  load  currents  and  radial  guide  surface  currents  are  given  in 
terms  of  E  and  B  fields  derived  from  the  transverse  axial  component  of  the 
generalized  Hertz  vector,  Z  ■  *5(r,t). 

The  next  step  in  a  complete  plasma/diode  model  is  the  constituitive 
relation  j_  -  OE;  The  Braginskii  transport  coefficients  express  currents  in 
the  plasma  load  in  terms  of  the  fields,  and  the  result  can  be  interpreted  as 
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requiring  (a)  radial  force  balance  on  an  "average  electron"  which  responds 

to  the  space-charge  field  E^  created  by  pressure  gradients  and  axial 

current,  and  (b)  constant  axial  electron  drift  speed  on  the  electron- 

inertia  timescale,  i.e.,  adiabatic  axial  electron  currents,  changing  on 

the  slower  timescale  of  the  fields.  Assuming  that  the  axial  drift  speed 

(Jz/ne)  and  the  radial  stress  quickly  relax  to  quasi-equilibrium  on  time- 

scales  fast  compared  with  those  for  tield  changes,  one  can  make  a  fully 

self-consistent  treatment  of  the  evolution  of  the  radial  fluid  velocity 

V  (r,t),  axial  current  density  J  (r,t)  and  the  incident  electromagnetic 
'  r  z 

field  2^(r,t)  simultaneously.  The  details  are  discussed  in  Appendix  C. 

The  theory  developed  above  is  an  adequate  basis  for  the  proper  calcula¬ 
tion  of  the  field/plasma  momentum  exchange  in  any  radial  implosion  whenever 
the  conductivity  relation  is  scalar  and  linear.  In  the  classical,  weakly 
coupled  plasma  of  course  this  relation  is  not  strictly  scalar,  but  to  a  good 
approximation  the  tensor  character  of  the  Ohms  law  can  be  superimposed  on 
the  electromagnetic  field  picture  developed  above  as  a  quasi-static  pertur¬ 
bation.  The  new  field  component  which  must  arise  (within  the  plasma  load) 
is  an  ambipolar  radial  E  field  due  to  radial  pressure  gradients.  Hall  effects, 
and  thermo-electric  effects.  In  the  present  case  this  radial  field  will  be 
independent  of  axial  coordinate  and  will  vanish  in  the  vacuum  region  of  the 
waveguide  because  it  is  essentially  the  field  of  a  radially  varying  axially 
uniform  line  image  charge  Qj(r)  with  Qj(r)  +  0  as  o(r,t)  +  0  in  r.  Using 
the  radial  stresses  developed  from  both  pE  I  ,  .  ,  and  aE  x  BQ  one  can 

therefore  define  a  completely  self-consistent  momentum  transfer  between  the 
far  waveguide  electromagnetic  field  and  the  annular  plasma  load.  The  result 
will  be  a  more  definitive  picture  of  the  load  dynamics  than  has  been  possible 
heretofore  in  the  classical  weak  coupled  plasma  limit.  This  detailed 
model  will  be  developed  further  in  Section  E,  but  a  simple  illustration 
of  the  physical  content  of  this  vector  Z(r,t)  will  serve  to  motivate  that 
analysis. 

C.  Simple  Solutions  in  the  Steady  State  Plane  Parallel  Waveguide 


In  the  case  of  steady  fields  th^  solution  for  Z(r,t)  assumes  a 
quite  symmetric  form  which  can  be  seen  to  be  a  firm  time  asymptotic  limit 
of  the  solutions  to  the  complete  wave  equation  B8.  Consider  a  constant 


conductivity  (oq)  ,  immovable  (Vp  -*■  0) ,  plug  (r  <  rp)  in  the  plane  parallel 
waveguide  and  transform  B8  using  the  dimensionless  variables 


r/rn  • 


ct/r 


o  ♦ 


Lm  r  a  /c 
o  o  ’ 


(12) 


based  on  an  arbitrary  scale  radius  r  ,  which  can  be  conveniently  set  as 
r  .  The  integrations  over  the  conducting  region  produce,  two  wave  equations 
with  continuous  sources  on  the  domains 


x  >  1:  x-1  3  x  3  Z  -  32Z 

xx  T 

x  <  1:  x-1  3  x  3  Z  -  32Z 

—  X  X  T 


and 


4*  JLnx  3t[Zq(x  3xz^J 

4t[£0  z(x>t>]“  \\£o  Zd.T)]) 


(13) 


The  steady  Ez>  Bg  fields  available  to  this  system  are  described  by  Z 
fields  which  arise  as  solutions  to  the  homogeneous  wave  equation.  The 
usual  vector  potential  Az(x,T),  cf  B9b,  exhibits  a  clear  decomposition 
into  Ez  terms  and  Bg  terms,  which  eliminates  the  scalar  potential  entirely 
as  previously  anticipated.  In  particular  one  has 

z(x,t)  *  Q0  (j-  +  t)  >  with 


(X>1) 

A*(x,x) 

*  r  ^  Q  (t+  2tt  X 
o  o  o 

x)  , 

Bg(x,x) 

211 <y 

cr  ’ 

(14) 

(x<l) 

Az(x,t) 

-  r;1  Qo  (T  +  TTZo(x2  - 

-  D) 

r 

2I(r<r  ) 

Bg(X.T) 

■  -  L2’  E=  VrJ  J  x  - 

_ E_ 
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and  Ez  =  -  QQ/r2  everywhere.  The  fields  here  corresponding  to  those  of 
a  leaky  capacitor  with  a  constant  voltage  source  at  x  -*•  00  maintaining  a 
steady  voltage. 

The  potential  clearly  decomposes  into  a  space-like  component ,  x2 , 
and  a  time-like  component,  T2.  The  space-like  function  represents  the 
charge  contained  at  any  time  within  the  radius  x  on  the  waveguide  plane, 
while  the  time-like  function  represents  the  cumulative  charge  involved 
in  the  continuous  superposition  of  incoming  and  outgoing  TEM  waves  which 
add  up  to  the  steady,  uniform  vector  field  Ez.  The  ordinary  vector 
potential  A also  exhibits  this  type  of  decomposition,  with  a  continuous 
value  and  spatial  derivative  at  the  conducting  interface  x  *  1 .  In  con¬ 
trast  to  a  more  conventional  formulation  of  the  same  problem  using  distinct 
scalar  and  vector  potentials,  the  present  method  makes  a  smooth  transition 
to  the  time  dependent  problem  posed  by  time  variations  in  either  Eq  or 

E  I  .In  these  situations  the  source  terms  are  forced  to  be  non-zero 

z'x  -*  °° 

because  3t(EqZ)  cannot  remain  spatially  independent  and  accommodate  the 

boundary  conditions.  Outside  rp(x  >  1),  the  time  dependent  value  of 

Eq(x  3^  Z|x  spawns  outward  running  waves  which  communicate  changes  in 
o 

the  conductor  to  the  charge  distributions  on  the  guide  planes.  Inside 
the  value  of  E  Z  at  the  interface  (x  =1)  acts  as  a  source  of  damped  waves 

o  -  p 

which  continue  to  penetrate  the  conductor  by  propagation  and  diffusion 
until  a  uniform  Ez  (free  wave)  solution  is  again  established  at  some  new 
boundary  value  Ezlx  ^  00.  Treating  the  same  or  a  similar  problem  using  a 
separate  set  of  A,  <J>  this  simple  process  of  charging  the  capacitor  or 
changing  the  conductivity  requires  some  connection  between  the  two  poten¬ 
tial  functions  in  order  that  Bg,  Ez  become  properly  interdependent.  Such 
a  connection,  or  the  complete  elimination  of  <f>,  is  not  always  transparent 
and  the  present  Hertz  vector  formalism  achieves  this  elimination  quite 
easily. 

It  might  appear  a  viable  conjecture  that  such  symmetry  breaking  as 
seen  in  the  time  dependent  "leaky  capacitor"  problem,  the  destruction  of 
the  space-like  plus  time-like  Hertz  potential  solutions,  is  a  fundamental 
property  of  the  theory.  This  is  in  fact  the  case,  as  discussed  in  the 
following  section. 


D.  Explicit  Source  Transformation  of  the  Wave  Equation  and  the 
Separation  of  Time  Scales 

The  homogeneous  relation  (B8)  which  specifies  the  evolution  of 
Z(r,t)  can  be  simplified  considerably  at  the  expense  of  introducing  a 
manifestly  physical  but  somewhat  complicated  source  term  involving  the 
current  density  profile  of  the  plasma.  Again,  setting  the  equation  in 
dimensionless  variables  and  choosing  an  (arbitrary)  dimensionless  scale 
radius  x^  in  the  static  Greens  function  (a  radius  sufficient  to  contain 
the  bulk  of  the  load  mass,  for  example),  the  terms  of  B8  which  contain 
all  of  the  source  information  exhibit  the  following  identity 


IZ  + 


OO  °°  /  ***  \  °o 

j  dx  Z-3^1  -  j  dx  3^.1  x£n^j  3^Z  -  Jln/^-j  '  f  dx  (S^O^z) 


-  j dx  x  Jin  E(x,t)  £°(x,t)  , 


(15) 


where  £°  =  -  x  H  x3  Z  and  x>  is  the  greater  of  x,  x.  When  this  result 

X  X  _ 

is  combined  with  the  c  v  (jc,t)  (=8r)  convection  terms  a  simple,  very 
physical  statement  obtains 

00 

3*Z  -  x_13x  (x3xZ)  =  4tt3t fdxx  Jin  Z(xx)  £  (3?,t)J  (16) 

o 

oo  ^ 

where  <£  =  £ °  -  8^3x  ("s^Z  -  J  dx  x  Hn  ^  Z  is  the  (dimensionless) 

L  o 

representation  of  the  field  E^  in  the  convecting  frame  Br.  The  generaliza¬ 
tion  to  all  orders  in  8r  is  quite  straightforward  (and  generates  the  space 
ordered  exponential  series  discussed  in  Appendix  B) ,  but  it  is  not 
-equired  when  the  convection  speed  implies  Sr  <<  1,  as  is  usually  the  case 
in  experimental  situations  of  interest.  This  new  relationship  displays 
quite  vividly  the  conjecture  of  the  previous  section.  The  source  term  is 
always  zero  in  any  situation  where  J  =  T.  £  is  time  independent  (for  any 
spatial  configuration  of  the  conductivity  Z)  and  thus  only  temporal  varia¬ 
tions  in  J  can  force  solutions  for  Z(x,x)  which  depart  from  the  simple 
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form  seen  in  the  previous  section  after  a  light  transit  time  across  the 
radial  domain  of  interest. 

A  second  result  of  (16)  above  is  that  a  clear  separation  of  the  time 
scales  occurs.  The  free  wave  solution,  being  characterized  by  the  light 
transit  time  across  the  system,  will  respond  to  the  much  slower  time 
variations  of  the  source  term  by  simply  propagating  variations  in  <^°  and 
J  over  the  system.  In  any  numerical  scheme  designed  to  solve  (16)  it  is 
therefore  essential  to  find  some  means  to  calculate  this  propagation  over 
large  intervals  in  t,  and  allow  a  concentration  of  the  numerical  effort 
on  tracking  the  temporal  variations  of  the  source.  In  Appendix  D  an 
explicit  forward  quadrature  for  the  free  wave  propagation  is  developed 
which  removes  this  problem.  The  method  is  related  to  the  Poisson  integral 
representation  of  two  dimensional  wave  solutions,  but  it  is  developed 
directly  from  the  free  wave  Greens  function  in  two  spatial  dimensions. 

The  result  of  each  forward  quadrature  is  the  specification  of  ^(XjTi+1) 

and  3  Z  (x.T  )  on  a  specified  (easily  rezonable)  mesh;  and  these 

i+1  3  11 

values  are,  of  course,  sufficient  to  specify  all  the  sources  required  to 

compute  the  next  forward  step  once  the  temporal  and  spatial  variation  of 

8  ,E  is  known.  In  the  fully  self-consistent  formulation  the  evolution  of 

Br,E  is  provided  by  a  similar  forward  advance  of  hydrodynamic  plasma 

variables,  but  the  method  developed  in  Appendix  D  allows  any  source  of 

these  fields.  The  development  of  the  Hertz  vector  sources,  on  the  right 

of  (16)  above,  can  be  cast  in  a  relatively  compact  set  of  operations  which 

proceed  from  Z.,  3  *Z ,  and  use  the  wave  equation  itself  to  eliminate 
J  T  J  A  ,  A 

higher  order  time  derivations  of  Z  when  they  appear  in  3T  etc. 

E.  The  1-D  Electromagnetohydrodynamic  Model  for  the  Plane  Parallel 
Waveguide 

In  the  numerical  method  a  more  natural  variable  is  a  *  x!  since 
this  removes  some  coefficients  in  the  differential  operators  which  are 
singular  at  x  -*■  0.  The  variables  a  and  x  are  used  interchangeably  in 
the  following  discussion;  and  Qq  is  a  reference  charge  value  which  allows 
Z  to  be  dimensionless  as  well.  Any  forward  advance  of  the  system  must 
begin  with  a  specification  of  *Zj  =  Z(XjT^)  and  ^Z^  =  3T  Z(x^T^)  together 
with  the  plasma  state  variables  * 


■’->r T*»  . 


••• 


the  ion  density,  ionization  state,  chemical  potential,  electron  and  ion 
temperature,  and  radial  convection  field. 

These  fluid  variables  and  their  gradients  provide  the  values  of 

i*  i. 

T  ,  n.  required  to  calculate  3  1  and  its  higher  time  derivatives;  the 
J  j  i  i»T 

Hertz  field  variables  Z.  and  Z.  provide  the  required  input  to  compute 
A  i  •  1  J 

3 T  <®  =  <£!  and  its  higher  time  derivatives.  For  example,  the  systematic 

Ti  J  % 

specification  of  the  time  derivatives  of  requires 


*  =  -  4  (3  a3  *Z)  ,  or  E  *  -  Q  r  ^  in  lab  coordinates, 

J  a  a  j  z  xo  o 


l€.  -  -  2a* 


\  [(3>j  -r 


iB,  »  -  2a*4 


[".^J  -  I-  ] 


(17) 


4(3  ad  iZ).  -  2a 4  i0,(32  a3  *Z) ,  -  *6  -23** 
J  a  a  j  j  a  a  j  j 


[“.^j  -  ^  A  ip«vi] 


and  these  in  turn  specify,  when  3  E  is  included,  the  complete  source 
term  for  the  next  time  advance  of  Z,  Z,  viz. 
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From  these  same  fields  the  fluid  stress,  which  provides  Bj ,  can  be  cal¬ 
culated  from  the  formulation  in  Appendix  C  and  the  implied  j  x  B  force 
density  — 

“  Vo5  ■ 

The  external  voltage  source  can  be  modeled  by  enforcing  a  boundary 

value  for  £  at  some  large  radius  (as  a  function  of  time)  and  the  current 

-  A 

drawn  by  the  model  plasma  is  simply  the  spatial  integral  of  E  <5  over  the 
domain.  Drift  speed  limits  on  current  conduction  can  be  imposed  quasi- 
statically  at  every  time  step  in  much  the  same  manner  as  one  develops 

the  Braginskii  conductivity  in  Appendix  C.  As  the  plasma  density  falls, 

A 

withc^  given,  there  will  always  exist  some  radius  for  which  the  drift 
speed  u  will  be  forced  to  exceed  the  sound  speed  c  ,  and  this  will 

Z  ^  A  S 

define  implicitly  a  (nonlinear)  conductivity  E((?)  which  satisfies  the 
~  a  £ 

constraint  jQE(£)'  e^egnIcg.  In  practice  this  limitation  may  drift 
in  or  out  in  x  as  the  £  profile  develops  and  it  will  generalize  the 
present  SPLAT  coronal  region  in  a  quite  physical  manner. 

The  calculation  of  the  usual  classical  heat  sources  and  the  radiative 
losses  can  be  folded  into  the  model  as  further  fluid  source  terms  in  any 
of  a  variety  of  well  known  methods  based  on  previous  work. 
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CHAPTER  IV.  IMPLICATIONS  OF  DISCHARGE  BEADING  WITH 
LOW-DENSITY  INTERCONNECTION  REGIONS 


A.  Introduction 


The  late-time,  fully-developed  stage  of  sausage-like  "beading"  of 
z-pinch  discharges  can  give  rise  to  enhanced  electron  heating  (and  thus 
enhanced  radiative  losses) ,  caused  by  Ohmic  anomalous  heating  in  the 
constricted  regions.  In  this  chapter  this  transfer  of  energy  from  conden¬ 
sations  of  magnetic  field  energy  to  radiation  is  examined  quantitatively, 
based  on  assumptions  about  the  nonlinear  state  of  the  instability. 


Constricted  portions  of  the  discharge,  with  low  density  and  cross- 
section,  and  mostly  or  entirely  "anomalous"  current,  alternate  with  the 
higher  density  "beads",  which  carry  current  classically.  The  extreme 
limit  of  this  phenomenon  is  that  of  multiple  diodes  in  series,  with  the 
nearly-evacuated  low  density  regions  considered  as  bipolar-flow  diodes, 
with  pinched  electron  flow,  as  suggested  by  Goldstein11.  In  all 
probability,  the  low-density  regions  cannot  evacuate  to  the  extent  required 
for  such  vacuum-diode  behavior. 

The  overall  resistive  heating  rate  is  of  course  VI,  with  the  total 
current  I  given  by  appropriate  circuit  equations,  but  the  local  heating 
rates  for  electrons  in  the  low  density  regions  are  balanced  by  increased 
radiative  loss  when  these  hotter  electrons  collide  with  the  denser  blobs 
of  plasma.  The  blobs  cannot  respond  hydrodynamically  to  the  increased 
heating  before  radiation  loses  the  deposited  energy. 
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B.  Current 


We  first  examine  the  conditions  under  which  all  the  current  in  the 

low-density  regions  would  be  drift-speed  limited,  i.e.,  (approximately) 

when  the  classical  drift  speed  =  oE/nee  would  exceed  the  sound  speed 

c  .  From  0 .  =  0)2/4i:v  .  and 
s  -*•  p  ei 

v^sec'1)  -  0.9  X  10"  (^j  (fj)’  (^t)  T-^2  (18) 

we  get 

•  1-8  »  10"  (^)  (f)’  E(MV/an)  (19) 

where  2  is  the  degree  of  ionization  (ng  =  zn^) ,  0j_  is  the  axial  conduc¬ 
tivity,  (across  magnetic  field  lines)  and  InA  is  the  plasma  parameter. 

The  classical  drift  speed  rises  rapidly  with  T  as  the  plasma  becomes  less 

e 

collisional. 


This  is  to  be  compared  with 


c 


s 


(20) 


For  10-times  ionized  aluminum  (A  -  26,  z  »  10),  one  sees  that  with 
electrons  at  keV  temperatures  the  classical  drift  speed  greatly  exceeds 
the  sound  speed  for  MV/cm  fields  unless  n^  is  of  order  1021  cm  or 
greater.  Elsewhere  we  argue  that  the  actual  drift  speed  cannot  much 
exceed  the  sound  speed,  so  that  the  current  is  drift-speed  limited  to  a 
value  approximately 


I  (MA) 
& 


(21) 


in  the  anomalous  ("a")  low-density  regions  of  cross-sectional  area  ffr2. 
Since  the  current  density  is  divergenceless,  this  must  be  the  current 
in  the  circuit. 
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C.  Voltage  Drop  in  Low  Density  Regions 


We  now  envision  a  number  N  of  such  millimeter  size  low-density 

"anomalous"  regions  of  length  alternating  with  denser  "classical" 

regions  ("beads")  of  length  ic.  Each  "classical"  region  has,  in  general, 

a  core  which  carries  no  current  because  it  has  not  been  penetrated  by  the 

magnetic  field  (and  hence  has  no  inductive  E) ,  a  classical  skin  current 

layer,  and  an  "anomalous"  (drift-speed  limited)  outer  corona.  The 

electrical  load  is  thus  envisioned  as  a  transmission  line  consisting  of 

N  elements  in  series,  with  each  element  as  in  Figure  12.  Since  the  total 

current  is  divergenceless,  it  is  given  by  Eq.(21),  and  depends  on  the 

applied  voltage  only  through  the  temperature.  If  the  corona  carries  only 

a  small  fraction  (1-8«1)  of  the  current,  the  DC  voltage  drop  along  each 

classical  bead  is  R  I  where  R  is  the  classical  resistance  of  the  bead, 
c  c 


i.e. 


R  (ft)-  3.5  x  10 

c 


£c(mm) 

ir(r£  -  rj)(mmz) 


(22) 


with  r  the  bead  radius  and  r . (<  r  )  the  current  penetration  radius. 

C  1C 

•  • 

One  can  see  that  when  inductive  (LI  and  LI)  effects  are  small  (i.e., 
if  L  *=  0  when  I  peaks  or  dips) ,  the  classical  resistance  of  a  few  tens  of 
1  mm  beads  at  roughly  1  keV  is  quite  small  if  the  current  penetrates  the 
beads,  and  almost  all  the  voltage  drop  occurs  in  the  anomalous  regions. 

On  the  other  hand,  r. /r  tends  to  be  frozen  in  at  early  times  by  the 

1  c  -1 
formation  of  the  high  conductivity,  and  tends  to  have  values  of  order  10 

or  less,  so  the  resistance  is  larger  than  the  current-penetrated  value  by 

an  order  of  magnitude  or  more.  Still,  this  gives  only  a  small  resistance 

for  the  "classical"  portioq  of  the  circuit  and  most  of  the  DC  voltage  drop 

would  occur  in  the  anomalous  regions. 
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D.  Joule  Heat ini 


Setting  aside  for  now  the  question  of  inductive  corrections,  this 
allows  us  to  estimate  the  Joule  heating  rate  in  the  anomalous  regions  in 
the  DC  case: 


E-I 


VI 

n  £ 

a  a 


N  £ 

N  £  Rc 
a  a 


(23) 


where  the  last  term  is  neglected  and  where  I  is  given  by  Eq.  (4) .  If  all 
this  energy  went  into  electron  temperature  and  none  escaped  to  ions  or 
photons,  it  would  correspond  to  a  rise 


•  2  VI 

T  a  —  — 

3  N  £ 
a  a 


irr‘  zn. 


i  •  e  4 , 


T (keV/ns)  at  15  V(MV) 


(24) 


(25) 


a  number  of  order  15  keV/ns. 


The  transit  time  of  a  typical  electron  from  one  bead  to  the  next  at 
a  constant  sound  speed  cg(^jtey)  would  be 

_q  /  a  10  -Is 

T  "  VCs  "  4,5  X  10  SeC  *  *aCnnn)  \26  T )  TkeV  (26) 

or  a  few  ns,  so  that  an  electron  coi'ld  acquire  several  tens  of  keV  before 
plowing  into  the  adjacent  dense  bead  and  losing  its  energy  to  dense- 
plasma  thermal  energy  and  not-so-cold-target  bremsstrahlung.  In  fact,  a 
certain  amount  of  this  energy  is  lost  to  ions  and  radiation  during  transit 
across  the  anomalous  region,  and  we  now  give  estimates  showing  that  loss  to 
be  dominated  by  classical  heat  conduction. 
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E*  Cooling  Rates  Offsetting  Joule  HeaHnfi 


The  classical  collisional  electron  cooling  time  to  ions  of  atomic 
mass  A  (taken  from  Book,  1980)  is 


T(keV/ns)  »-!■(.  10'2  (fj)’  (*§*)  (?)  (j^f) 

which  is  negligible  at  1  keV  and  above  (unless  ng  >  102 1  cmJ ) . 
radiative  cooling  gives  approximately 


T(keV/ns) 


(27) 


The 


(28) 


for  aluminum  with  1  keV  (Terry  and  Guillory,  Ref.  1  p.  5  and  plasma 

radiation  transport  computations  done  by  Duston  (1980)). 


The  classical  heat  flux,  which  is  dominated  by  the  Hall  flux  term 

^  x  V.T  since  v  /u >  is  small  except  very  near  the  axis,  would  give  a 

jL  e  e  ce 

cooling  rate 


T(keV/ns) -  3.5  x  10 


r 

a 


(29) 


where  r^,  is  the  radial  temperature  scaleheight. 


12  13 

But  the  modified  two-stream  instability  ’  grows  at  a  rate  about 


half  the  lower-hybrid  frequency 


.3  x  10 


1 1 


/26\  (zj  (  ni  )(  ra  V* 

\A  /  \10/  \10l •/ \1  mm  / 


lkeV 


(30) 


(based  on  Eq.  (21)  for  I)  and  thus  can  easily  saturate  in  less  than  1  ns, 
leading  to  alteration  of  the  diffusive  heat  flux  as  well  as  being  the 
source  of  the  drift-speed  limiting  assumed  at  the  outset. 
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But  the  alteration  of  the  heat  flux,  according  to  the  form  of 

14 

anomalous  heat  conductivity  given  by  Manheimer 


•cT/16  eB 


(31) 


is  small  and  gives  a  contribution 


T(keV/ns)  ~  10 


-3 


/AN5*  /10\ 

\26/  U  ) 


3/2 


/1018\  f  r(mm)  1 

\  n±  )  [ffr|(mm2K2  (mm2)  J 


.-3/2 

keV 


(32) 


to  the  electron  cooling  in  transit  across  the  anomalous  low-density  regions. 
(The  inverse  dependence  on  n^  comes  only  from  the  linear  dependence  of  I, 
and  thus  B,  on  n^.) 

The  instability  tries  to  stabilize  by  bringing  a  small  number  of 
ions  up  to  the  sound  speed^  but  this  is  an  important  energy  drain  for  the 
electrons  only  where  T  »  T . .  For  the  time  development  of  T  due  to  this 
anomalous  heating,  one  may  consult  Lampe  et  al,  .  A  fraction  fi~(me/mi) 
of  the  ions  is  accelerated  to  a  speed  comparable  with  c  ,  i.e.,  to  energies 

3  3 

-j  Te>  (The  directions  are  isotropized  by  the  magnetic  field.)  There  is 
thus  a  loss  of  ion  thermal  energy  due  to  the  exit  of  these  particles  from 
the  low  density  region,  which  is  only  partly  compensated  by  the  thermal 
influx  of  fons  at  temperature  T^.  The  net  heat  loss  flux  per  cm2  of  area  is 


"ifi[(K)(K)-(ivThi)  (t  ti)]  ( 

(vThi  is  the  ion  thermal  speed);  and  when  the  instability  is  saturated 
and  quasi-steady  this  energy  loss  is  made  up  by  energy  transfer  from 
electrons  to  ions  via  unstable  waves: 

|  Ve  *ri  *a  "  *rl  (I  nifi)[csTe  -  vThi  Ti]  • 
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This  gives  rise  to  an  electron  cooling  rate 


*.<«/->  i  -  *  “-3  m"  {?r  • 

This  rate  is  small  because  of  the  small  number  of  ions  involved, 
relative  to  the  number  of  electrons.  In  fact,  for  nominal  parameters 
all  of  the  electron  cooling  rates  (27)  -  (34)  are  small  compared  with 
the  Ohmic  heating  rate,  and  so  most  of  the  Ohmic  energy  gain  is  delivered 
to  the  blobs  after  transit  across  the  anomalous  region. 
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F .  Power  Deliver' 


Assuming  the  dense  blobs  to  be  thicker  than  one  range  for  the  15-50 
keV  electrons  accelerated  into  them,  the  total  net  power  delivery  to 
them  is 


(T  -  T  )  NUT2 
2  e  eo  a 


(35) 


where  Tg  -  Tgo  is  the  temperature  gain  by  Ohmic  heating,  less  losses, 
during  transit  from  blob  to  blob.  If  the  cooling  during  transit  is 
negligible  and  if  almost  all  the  voltage  drop  occurs  across  the  low-density 
regions,  then  the  power  delivered  is 


P(T“>  "  °'36  (l5»)  fe)  (lo)3/2  (ff  (w>)  4 


with  n^,  representing  either  their  average  values  or  their  initial 

(unheated)  values  in  the  low-density  region.  If  less  than  the  full  voltage 
drop  occurs  across  the  anomalous  regions,  the  value  of  V  above  is  adjusted 
downward  accordingly. 


G.  Peak  Electron  Temperature 

In  the  low-density  regions,  the  electron  temperature  varies 
axially,  increasing  in  the  direction  of  electron  flow  because  of 
collisional  and  unstable-wave  isotropization  of  the  energy  gained 
from  the  axial  electric  field.  The  temperature  begins  at  approxi¬ 
mately  its  value  in  the  adjacent  classical  bead  (a  quasi-bala;ice  of 
deposition  and  radiative  loss),  increases  as  the  electron  traverses 
the  low  density  region,  and  reaches  its  maximum  just  before  deposition 
in  the  next  "bead"  downstream. 

•  i< 

Since  the  heating  during  transit  has  the  approximate  form  Te® 
(from  Eq.  (25),  one  has  for  the  case  of  negligible  heat  conduction  and 
radiation 


cg(t) ^  cg(0)  +  cgt 
z(t)  =  Cg (t) 


(37) 


So 


that  z(t)  =  c  (0)t  +  I  c  tz  and  the  final  temperature  is  such  that 
s  2  s 


c  (t.)  =  Sc1  +  2  c  £ 
s  f  so  s  a 


i.e  T  (£  )  -  T  + 
e  a  eo 


670  keV 
N 


jv(MV)  -  Llj 


(38) 

(39) 


from  Eqs.  (20)  and (25),  as  an  overestimate  of  the  peak  temperature. 

(N  ,  the  number  of  anomalous  regions,  is  observed  to  be  typically 
of  the  order  of  20  to  30,  while  the  inductively  corrected  voltage 
V  -  LI  is  one  or  two  MV.)  The  voltage  drop  along  the  classical  beads 
has  been  neglected  compared  with  that  over  the  anomalous  regions 
between  beads,  and  including  it  reduces  the  effective  voltage  appearing 
in  Eq.  (39)  to  that  occurring  only  over  the  Ng  anomalous  regions.  In 
this  approximation  the  temperature  increases  linearly  with  distance 
across  each  acceleration  region. 
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The  indication  that  the  peak  electron  temperature  may  be  tens 
of  keV  implies  that  Bremsstrahlung  and  inner  shell  x-rays  can  be  produced, 
predominantly  on  the  cathode  ends  of  the  classical  beads  if  their  density 
is  high  enough,  or  throughout  the  beads  if  they  are  only  about  one  range 
thick  for  the  energetic  electrons. 
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r^X-riCMi*k' 


H.  Inductive  Effects 


The  generalized  Ohm's  lav  in  a  classical  medium  takes  the  form 


E«nJ  +  kj--xB 
e 


or  V  -  RI  +  LI  +  LI  .  (41) 

(k  is  a  constant  related  to  the  electron  inertia  and  the  geometry,  and 
L  is  the  self-inductance.)  Voltages  induced  by  changing  the  current- 
carrying  radius  r&  are  of  order 

*a  -1 

E  (MV/cm)  ^  21  (MA)  ^  (ns  A)  ,  (42  } 

a 

which  becomes  of  order  1  MV  just  before  the  peak  compression.  The 
corresponding  impedance 

L  (Ohms)  ~2£(cm)  r  /r  (ns  *)  (43) 

a  a 

has  maximum  value  on  the  order  of  1  Ohm,  in  series  with  1-4  Ohms  effec¬ 
tive  resistance,  and  so  is  not  a  negligible  correction  except  for  a  few- 

nanosecond  intervals  when  r  0.  In  SPLAT  code  runs,  the  L  term  is 

a 

even  somewhat  larger,  at  its  maximum,  than  the  Ohmic  resistance  of  the 
classical  plasma  core .  But  where  the  current  becomes  limited  by  anomalous 
regions  in  series  with  classical  ones,  the  effective  resistance  R*  in 
the  equation 

V  -  R*I  +  LI  +  Li  (44) 

is  larger  (and  voltage-dependent,  as  V);  so  the  L  term  no  longer 
dominates,  although  it  is  not  strictly  negligible  except  very  near  peak 
compression  and  for  the  "pause"  phase  of  "pause-and-collapse"  compressions. 
For  our  simple  description  we  neglect  the  L  term,  which  is  acceptable 
when  V/I  is  much  larger  than  L,  i.e.. 


4 

The  LI  term  may  be  governed  by  nonlinear  sausage  instability 

dynamics  not  yet  well  understood;  when  I  is  given  by  Eq.(21),  I  is 

proportional  to  4— (n.nr2) ,  i.e. ,  the  change  in  particle  number  in  the 
at  i  a 

necked-off  region.  This  subject  is  to  be  investigated  uner  1981  and 
1982  funding,  and  is  beyond  the  scope  of  available  theory  at  the  time 
of  this  report. 


I .  Justification  of  Drift-Speed— Limited  Current 


When  the  ion  and  electron  temperatures  are  comparable  one  does 
not  get  appreciable  growth  of  the  ion  acoustic  instability  (Kovrizhnykh, 
1967;  Biskamp  and  Chodura,  1971,^),  but  the  modified  two-stream  insta¬ 
bility  (Ott  et  al;^  McBride  et  al,2^)  grows  with  linear  growth  rate 


Im  u)(MTS) 


TU.,.  =  -W  U  .(1  +  0)2  /tO2  )~^2 

2  LH  2  pi  pe  ce 


1 

2 


(46) 


where  o>TH  is  the  "lower  hybrid"  frequency,  ojpe  is  the  plasma  frequency, 

and  u>ce  is  the  electron  gyrofrequency .  This  persists  as  long  as  the 

drift  velocity  v^  exceeds  approximately  the  sound  speed,  cg,  and  is 

marginally  stable  when  v.%c  .  The  nonlinear  saturation  time,  of 
-1  as 

order  300  u>  ,  is  much  shorter  than  the  £/c„  transit  time  of  electrons 
pe  s 

across  the  anode-cathode  gap,  so  the  waves  can  grow  to  saturation  within 
the  diode. 


Including  the  electron  magnetization  more  carefully,  one  also  gets 

21 

an  "electron  cyclotron  drift  instability"  (Forslund,  Morse  and  Nielsen 
with  growth  rate 


Im  oi(ECD) 


ce  v 


e 


ei 


(47> 


where  v^  is  the  electron  drift  speed,  vg  is  the  electron  thermal  speed, 
and  the  electron  collision  frequency  as  before.  But  for  typical 
parameter  values,  v<j/vg  should  be  of  order  /A/Z  so  the  two  linear 
growth  rates  are  comparable: 


Im  u 


2.6  x  1011 


I  (MA) 

r  (mm) 
a 


as  long  as  this  exceeds  in  Eq.  (18). 


(48) 
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The  radial  gradients  in  magnetic  fields  can  also  give  rise  to  a 
"lower-hybrid  drift"  instability  (Krall  and  Liewer,^)  with  linear 
growth  rate 


Im  w(LHD) 


~  “pi[5B 


v2  1  % 
dB  ___e_  I 

dr  u)  cl 
ce  s J 


(id  /id  ) 
ce  pe 


7.4  x  10 


[?  (r)' 


h2B(MG)2 


where  h„ 


=  (b  df)  ~ra(an)* 


Since  B  is  at  most  2I/cra>  this  gives 

I"  “(U,D)  •<-  10”  ^=o[(f)(f )’  U1M  ($ 


which  is  generally  less  than  the  MTS  and  ECD  instability  growth  rates. 

Because  the  modified  two-stream  instability  exponentiates  in  a 
few  picoseconds  and  saturates  in  several  growth  times,  it  is  most  likely 
that  the  electron  drift  speed  is  limited  in  this  way  to  values  near 
that  which  stabilizes  the  instability,  namely  the  sound  speed.  This 
is  the  justification  for  the  anomalous  drift-speed-limited  current- 
carrying  behavior  assumed  earlier. 
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J.  Application  to  Cylindrical  Implosions  and  the  “SPLAT” 


From  Eqs.  (19)  and  (20) ,  we  easily  derive  a  fact  of  some  significance. 

Before  the  field  has  penetrated  a  cylindrical  z-pinch  plasma,  the  radial 

profile  of  the  E  electric  field  of  course  goes  from  negligible  values  at 
z 

the  current  penetrsftion  front  to  values  of  order  V/fc  at  large  radii.  For 
all  radii  where 


E(r,  (Mv/„)  >  ^  (f§4)fe; 


,5/2 


) 


(50) 


the  current  is  drift-speed  limited.  For  temperatures  of  order  keV  and 
above,  the  critical  E(r)  value  at  r*a  is  much  less  than  V/Jt  'v*  lMV/cm  when 


n^(r)  is  small  compared  with  102°. 


>-7 


But  from  the  voltage  rise  time  Ty  10  sec,  we  can  estimate  a 


classical  skin  depth  in  the  dense  (classical)  region: 


.-3/4 


A  (cm)  -  29  /rv(sec)  TkeV 


/ £nA  z  \ 

^10  10 ) 


(51) 


and  this  gives  a  core  current, 
a 


-(a-r) /A 


I  J  2nraE(a)e  *"“dr,  i.e., 
o 

A 


[j(MA)  -  0.21  y!  T^4  y)  (f)  .<«.)  ^ 


n<(a) 
n~ 


(52) 


Uhen  the  core  carries  most  of  the  current  (1^*1^) ,  this  gives  for  the 
density  at  the  critical  surface  of  radius  a 


-  4  R  ^-15-  aWioVt*  ilMl 

Qi»  4,8  \ln\  1GJ  (  z  J  TkeV  y  xy  j  a(nin) 


10 


(53) 
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Thus,  for  a  keV  aluminum  plasma  carrying  1  MA  current  at  the  skin  of  a 
1  mm2  area,  if  the  on-axis  density  drops  below  about  1019  ions/cm3  any¬ 
where,  then  the  classical  "core"  tends  to  disappear  there.  And  since 

2  0 

this  value  (Eq.  53)  of  n^  tends  to  be  small  compared  with  10  cm  ,  most 
of  the  E  field  is  shielded  from  the  classical  region,  i.e.,  E(a)«V/£, 
while  there  is  a  classical  region. 

The  one-dimensional  strongly-coupled  plasma  implosion  code  SPLAT 
(Terry  and  Guillory, *) sets  (Initially)  a  radius  a  beyond  which 
the  current  is  drift  speed  limited  at  the  initial  temperature.  In  the 
present  form  of  the  code,  this  radius  is  then  followed  in  a  Lagrangian 
manner,  moving  with  the  fluid  velocity.  In  fact,  however,  as  the 
temperature  and  current  increase,  the  real  "critical  surface"  at  which 
conduction  becomes  anamalous  moves  inward  more  rapidly  than  the  particles, 
following  the  variable  density  contour  of  Eq.  (53).  Thus  a  larger  and 
increasing  fractionof  the  current,  Ij/I*  is  carried  "anomalously."  This 
fraction  is  thus  underestimated  by  the  code.  The  corona  is  actually 
more  important  than  indicated  by  the  code  for  another  reason:  its  density 
at  the  true  "critical  surface"  r=a  is  larger  than  presently  accounted 
for  and  the  total  ohmic  heating  should  be  larger. 


K.  Sausage  Instability  Dynamics 

The  nonlinear  behavior  of  the  axisymmetric  sausage  instability 
has  been  treated  only  in  a  few  limiting  cases.  A  nonradiating, 
perfectly  conducting  incompressible  MHD  model  with  uniform  density 
and  no  heat  conduction  was  treated  by  Book,  Ott  and  Lampe.1^ 

In  both  the  long  wavelength  (X  >>  r)  and  short  wavelength  limits, 
these  assumptions  gave  tractable  equations,  which  showed  a  spool-like 
development  of  the  plasma  shape,  with  thin  rapidly  expanding  disks 
of  enlarged  radius  and  broad,  slowly  constricting  regions  of  reduced 
radius.  The  presence  of  finite  conductivity  and  axial  electric  field 
certainly  modifies  the  dynamics;  and  even  more  so,  the  inclusion  of 
finite  compressibility,  which  is  being  studied  now  by  Book  et  al,  may 
greatly  alter  the  conclusions  made  in  the  simpler  model. 

Because  of  this,  we  have  proposed  in  our  10/81-10/82  program  an 
extension  of  the  Book  et  al  theory  to  include  finite  conductivity  and 
E  field,  and  the  construction  of  a  2-dimensional  r-z  code  for  modeling 
the  behavior  with  radiation,  circuit  model,  etc.  coupled  into  the 
dynamics . 

The  equations  describing  the  motion  are: 


continuity: 

9tP  +  9z<pvz>  +  r  9r^PrvP  *  0 

pDv/Dt  -  VP: 

pOtvz  +  vz9zvz  +  vr3rvz)  -  -  9zP 

p(3tvr  +  vz9zvr  +  vr3rvr)  *  -  3rp 

state: 

P  -  yTp/mi  +  B2/81T  +  KB2/r  B2  -  B2(r,z) 

VxB: 

79r(rB)  -  —  V  -  3ZB  -  —  Jr 

VxE: 

9zEr  -  9rEz  “  -  i  9tB 

V*J«0: 

h  (rJ  )  +  3  J  -  0 

y  j  '  j-'  2  Z 
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Heat  balance:  T  specified  as  function  of  discharge  voltage 
and  p,  in  simplest  model,  or  governed  by  time 
dependent  heating  and  cooling  rates,  in  more 
complex  model. 

The  radial  component  of  J  cannot  be  neglected  because  then  3  B-0 

z 

does  not  allow  the  instability.  The  constants  k^  in  the  unmagnetized 
conductivity  and  kg  in  the  Hall  conductivity  can  be  found  in  the 
standard  textbooks;  the  constant  kj.  (see  Eq.  (21))  gives  the  drift-speed 
limited  current.  The  constant  K  provides  the  magnetic  hoop  stress  due 
to  3-line  curvature.  The  mass  density  is  denoted  by  o,  electron  (and 
ion)  temperature  by  T,  and  fluid  velocity  by  (v  ,v  ).  All  field  and 
fluid  quantities  are  functions  of  r  and  z.  Imposing  self-similar 
convex  profiles  of  p  with  radius  (e.g.,  Bennett  profiles)  and 
periodicity  in  z,  with  only  the  on-axis  density  and  Bennett  radius 
varying  in  z,  may  allow  certain  artificial  but  useful  simplifications. 
Making  the  ideal-MHD  assumption  is  not  admissible  here,  as  it  allows 
only  surface  currents  and  cannot  predict  the  density  drop  in  the 
necked-off  regions.  The  fluid  model  described  may  break  down  when 
these  regions  produce  non-stagnated  electron  flows. 
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CHAPTER  V.  CURRENT  PROBLEMS 


A.  Improved  numerical  methods 

As  pointed  out  above,  there  are  several  advantages  to  the  use  of  the 

generalized  Hertz  vector  Z(r,t)  as  a  potential  for  the  waveguide  fields 

driving  an  imploding  plasma.  The  price  of  this  smooth  generality  is  that 

all  the  observable  fields  are  two  derivatives  removed  from  the  potential. 

Hence  we  require  very  accurate  numerical  differentiation  methods  for  the 

estimation  of  Bg  from  any  discrete  representation,  i.e.  Z^i^)  and 

ZCXjT^).  Moreover,  as  the  plasma  load  absorbs  and  reflects  the  incoming 

electromagnetic  wave,  the  rapid  decay  of  E  (x,x)  within  the  plasma  demands 

z 

the  resolution  of  the  derivatives  of  Z(x,t)  over  a  very  small  spatial 
scale.  Adding  to  the  challenge,  one  must  assess  third  and  fourth  order 
derivatives  in  calculating  the  source  terms  in  the  wave  equation  for 
Z (x ,t) , 

The  usual  finite  difference  methods,  based  on  polynomial  interpolation, 

tend  to  produce  rather  noisy  estimates  of  higher  derivatives  unless  the 

22 

spatial  mesh  is  quite  dense.  Recent  work  by  Talmi  and  Gilat  ,  however, 
offers  a  solution  to  this  problem  by  using  a  so-called  smooth  interpolation. 
In  classical  analysis  the  error  functional,  measuring  the  global  deviation 
of  an  interpolant  from  the  original  function,  is  usually  based  on  the  Hilbert 
norm.  Minimizing  error  with  respect  to  this  norm  does  not  constrain  the 
derivatives  of  the  interpolant  and  often  leads  to  rather  bizarre  behavior. 

The  generalization  due  to  Talmi  and  Gilat  which  leads  to  a  smooth  Inter¬ 
polant  is  the  invention  of  a  norm  functional  which  involves  all  higher 
derivatives  as  well  as  function  amplitudes  in  computing  an  error  measure. 

This  new  method  has  been  implemented  recently  and  has  generated 
encouraging  results  interpolating  simple  test  functions.  Function  values 
can  be  represented  to  7  significant  figures,  and  second  derivatives  to  4 
significant  figures.  This  sort  of  performance  is  not  quite  sufficient  for 
use  with  the  Hertz  potential,  but  it  is  quite  adequate  for  the  representa¬ 
tions  of  CRE  radiation  and  ionization  dynamics  required  in  section  B,  below. 
The  Investigation  at  present  is  focused  on  the  optimal  tradeoff  between  the 
number  of  sample  points  and  the  correlation  width  which  determines  the 
effective  coupling  domain  among  points  on  the  mesh.  There  is  every  reason 
to  believe  that  the  satisfactory  representation  of  third  and  higher  deriva¬ 
tives  is  within  the  capabilities  of  this  new  method. 


B.  Alternatve  thermodynamic  variables  for  ionization  dynamics 

In  hydrodynamic  calculations  requiring  ionization  dynamics  a  central 

problem  is  the  branching  ratio,  for  external  energy  input  to  the  fluid, 

between  electron  heating  and  further  ionization.  As  each  new  atomic 

shell  structure  is  opened  to  ionization  (because  the  electron  thermal  energies 

and  radiation  intensity  have  achieved  an  appropriate  threshold)  this 

branching  ratio  changes  radically  with  a  large  fraction  of  the  input  energy 

going  to  further  ionization  rather  than  electron  heating.  On  the  other 

DT 

hand,  any  attempt  to  calculate  -^e  and  advance  the  plasma  state  with  ^ 

and  eT  (the  chemical  potential)  held  fixed,  only  to  then  require  a  fresh 

value  of  and  ^  at  some  new  temperature  Tg,  will  lead  to  serious  inaccuracy 

and  perhaps  numerical  instabilities  if  the  ionization  energy  is  recovered 
* 

from  Tg.  This  is  why  the  branching  ratios  3^  <1 ,  3^  and  (to  a  lesser 

,  .  \  ^  «\  _  j  .  .  6  6 


extent)  3  2,  3  £_  are  so  important. 

nI  Q  nx  I 


The  calculation  of  these  ratios  from  first  principles  within  the  CRE 
framework  is  not  an  easy  task,  but  their  estimation  by  means  of  the  smooth 
interpolation  discussed  above  is  both  simple  and  accurate.  The  CRE  model 
implies  an  effective  equation  of  state  £j(Te,  n^)  and^  (T^nj)  which  can  be 
used  to  simplify  the  simultaneous  calculations  of  electron  heating  and 
ionization.  Instead  of  using  Tg  one  may  define  a  new  energy  variable 


l0  =  +  £t(T  )  and  T  =  T  (9  ), 

2  e  29  e  I  e  e  ee 


[I,  il+  3  , 

d£T  ' 

2  »-  — J I 

+  2aj  .[It  |2  +^L  1 

2  eW  2  < 

L  e 

J  3T 

e  - 

Dt  1.2  e  Sn^  dn^  J 

its  inverse  function.  The  appropriate  hydrodynamic  equations  now  advance 
0g  and  Tj.  using  the  usual  heat  sources  and  i luxes;  viz. 


3  D_  e  _  DT 
2  Dt  e  Dt* 


^•X)-!-a-:IY  +  ^2(Te(ee)-TI).TIe  ( 

where  the  partial  derivatives  and  functions  9  (T  )  arising  from  the  CRE 

e  e 

model  are  represented  by  a  smooth  interpolant. 

Present  efforts  in  achieving  this  reformulation  are  directed  toward 
finding  the  best  temperature  and  density  domain  grid  and  optimal  inter¬ 
polation  parameters.  In  addition,  the  restructuring  of  present  hydro- 
dynamic  codes  to  accept  these  new  variables  is  underway. 

62 


REFERENCES 


1.  R.  E.  Terry  and  J.  Guillory  JAYCOR  Report  No.  TPD200-80-001-DFR, 

July  1980. 

2.  D.  Duston  and  J.  Davis,  Phys.  Rev.  A  21_,  May  1980. 

3.  D.  Duston,  private  communication. 

4.  J.  Apruzese,  private  communication. 

5.  R.  E.  Terry  and  J.  Guillory  JAYCOR  Report  No.  TPD-200-80-001-DFR, 
pp.  67-70,  July  1980. 

6.  M.  Gersten,  private  communication. 

7.  J.  Katzenstein,  "Optimum  Coupling  of  Imploding  Loads  to  Pulse 
Generators."  ,  Maxwell  Laboratory  Report. 

8.  R.  E.  Terry  and  J.  Guillory,  JAYCOR  Rpt.  No.  TPD200-80-012,  10/80,  pp.  13-15. 

9.  J.  A.  Stratton,  Electromagnetic  Theory,  McGraw  Hill,  Hew  York,  1941 
1st  Ed. ,  p.  28. 

10.  W.  R.  Smythe,  Static  and  Dynamic  Electricity,  McGraw  Hill,  New  York 
1968,  3rd  Ed.,  p.  416. 

11.  Shyke  Goldstein,  in  H.  W.  Bloomberg,  et  al.,  "Energy  Coupling  Mechanisms 
in  Multiple  Wire  Array  Configurations."  Science  Applications,  Inc. 

Report  #SAI-102-72-004 ,  (1979). 

12.  N.  A.  Krall  and  P.  C.  Liewer,  "Low-Frequency  Instabilities  in 
Magnetic  Pulses,"  Phys.  Rev.  A4,  2094  (1971). 

13.  E.  Ott,  J.  B.  McBride,  J.  H.  Orens,  and  J.  P.  Boris,  "Turbulent 
Heating  in  Computer  Simulations  of  the  Modified  Plasma  Two-Stream 
Instability,"  Phys.  Rev.  Lett.  28^  88  (1972). 

14.  W.  M.  Manheimer,  "Anomalous  Transport  from  Plasma  Waves,"  Journal  de 
Physique,  4£,  C7-269  (1979)  Suppl. ; 

15.  G.  E.  Vekshtein,  D.  D.  Ryutov  and  R.  Z.  Sagdeev,  "Asymptotic  Solution 
of  the  Problem  of  the  Anomalous  Resistance  of  a  Collisionless  Plasma." 

Sov.  Phys.  JETP  33,  1152  (1971). 


63 


REFERENCES  (cont.) 


16.  M.  Lampe,  W.  M.  Manheimer,  J.  B.  McBride,  J.  H.  Orens,  R.  Shanny 
and  R.  N.  Sudan,  "Nonlinear  Development  of  the  Beam-Cyclotron 
Instability,"  Phys.  Rev.  Lett.  26 1221  (1971). 

17.  D.  L.  Book,  E.  Ott,  and  M.  Lampe,  "Nonlinear  Evolution  of  the  Sausage 
Instability"  Phys.  Fluids  19,  1982  (1976). 

18.  L.  M.  Kovrizhnykh,  "Nonlinear  Theory  of  Current  Instability  in  a 
Non-isothermal  Plasma."  Sov.  Phys.  JETP  2k_,  1210  (1967). 

19.  D.  Biskamp  and  R.  Chodura,  "Computer  .Simulation  of  Anomalous 
Resistance."  Phys.  Rev.  Lett.  27_,  1553  (1971). 

20.  J.  McBride,  E.  Ott,  J.  Boris,  and  J.  Orens,  "Theory  and  Simulation 
of  Turbulent  Heating  by  the  Modified  Two-Stream  Instability."  Phys. 
Fluids  15,  2367  (1972). 

21.  D.  W.  Forslund,  R.  L.  Morse,  and  C.  W.  Nielson,  "Nonlinear  Electron- 
Cyclotron  Drift  Instability  and  Turbulence"  Phys.  Rev.  Lett.  27 , 

1424  (1971). 

22.  A.  Talmi  and  G.  Gilat,  Journal  of  Computational  Physics  23^,  93, 
(1977). 

23.  S.  I.  Braginskii  in  Reviews  of  Plasma  Physics,  ed.  M.  A.  Leontovich 
(Consultants  Bureau,  N.Y. ,  1965),  p.  249  ff. 

24.  P.  J.  Davis  and  I.  Polonsky,  Handbook  of  Mathematical  Functions, 

NBS  AMS-55  (ed.  M.  Abramovitz  and  I.  A.  Stegun) ,  USGPO  (1964). 

25.  A.  H.  Stroud,  Approximate  Calculation  of  Multiple  Integrals,  Prentice 
Hall,  Englewood  Cliff,  NJ  (1971),  (pp.  267-289). 


64 


APPENDIX  A 


FIELD  REPRESENTATIONS 

The  procedure  for  choosing  Z_  can  be  formalized  easily  and  points 
out  the  role  of  a  Hertz  vector  potential  in  the  conventional  scheme  of 
electrodynamics . 

Field  Representation  Theorem: 

a)  Given  potentials  A,  <p  satisfying  the  Lorentz  gauge  condition, 

d  c 

the  choice  Z^  _  cl  dt  A(jc,  t  )  implies  that  $  =  -  V*Z  and  (obviously) 

-1  ~  J  1  1 
that  A  =  c  3^. 

b)  Given  the  vector  field  Z,  the  choice  A  ^  c  ^  3^Z  and  <p  =  “  Z'Z 
implies  that  A,  <p  satisfy  a  Lorentz  gauge  condition.  Moreover, 
the  transformation  Z  -*■  Z'  =  Z_  +  with  p  any  solution  of  Q]2^=0  is 
the  generator  of  all  restricted  gauge  transformations  (within  the 
Lorentz  gauge)  which  leave  the  observable  fields  E^,  B  invariant. 

This  result  is  readily  established,  as  summarized  below.  With  A,  4>  in  the 
Lorentz  gauge  and  _Z  defined  through  A,  the  gauge  condition  on  A,<|>  implies 
that  4>  =  -  V’Z  and,  of  course,  A  =  c  ^  S^Z.  On  the  other  hand,  given  Z^ 
and  A,  <$>  defined  through  Z,  the  Lorentz  gauge  condition  is  an  immediate 
consequence,  and  simple  substitution  verifies  that  Z'  as  given  preserves 
the  gauge  and  leaves  E^,  j5  invariant.  Of  course,  if  A,  <j>  are  to  be  repre¬ 
sented  on  all  spacetime  by  one  field,  Z_,  it  must  be  shown  that  only  one 
non-trivial  propagation  equation  is  needed  for  evolving  this  new  field  Z_. 
This  is  in  fact  the  case. 

Uniqueness  Theorem  for  Z(x,t):  Let  initial  fields  E°  and  B°  be 
specified  on  all  space  at  some  time  tQ,  together  with  a  source  field 
T  (equivalently  p  and  J)  on  a  domain  (x,  t  tQ)  and  boundary  condi¬ 
tions  of  the  form  c^  E(x,t)  +  j3(x,t)  on  any  surfaces  of  discontin¬ 

uity  that  are  present.  The  generalized  Hertz  vector  Z(x,t),  as  a 
solution  of  V(V- Z)  -  VxVxZ  -  c  “S2  Z  =  -  4ttT,  is  then  specified 


A1 


uniquely  (to  within  a  trivial  restricted  gauge  transformation) 
and  stably  on  (x,  t  >  t  ). 

—  O 

The  proof  is  somewhat  lengthy  but,  in  summary,  decompose  Z  and  T  into 
C  t  £  ff 

solenoidal  (Z  ,  T  )  and  lamellar  (.Z  ,  T  )  parts.  For  the  solenoidal  fields, 

the  initial  specification  of  Et,  B  ,  and  JC(=  is  equivalent  to  an 

— o  — o  —  t — 

initial  specification  of  ±\t  and  3fcgt|  .  For  the  lamellar  fields  the 

o  o 

a  z  z 

initial  specification  of  E  (=  -  4ttT  )  and  J  is  equivalent  to  an  initial 
Z  °  o  ° 

specification  of  Z  \  and  3t  ZT  \ t  .  The  homogeneous  solution  for  must 

co  o 

then  be  selected  to  produce  the  specified  boundary  conditions.  The  speci¬ 
fication  of  zj  and  3^)  required  by  the  initial  conditions  thus  consti- 
co  o 

tutes  Cauchy  open  surface  boundary  conditions  and  a  unique,  stable  solution 
exists  for  the  Z_  wave  equation.  The  ambiguity  in  Z  with  respect  to 
the  restricted  gauge  transformations  can  be  erased  by  simply  setting  to 
zero  all  trivial  constants  of  integration  in  the  inverse  mapping  (E°,  B°)  ■* 

(Z|t  .  3tilt  }• 

o  o 


A2 


APPENDIX  B 


HERTZ  VECTOR  POTENTIAL  FOR  THE  PLASMA  LOADED  WAVEGUIDE 


Using  the  Green's  function  one  can  decompose  any  vector  field  b(x,t) 
into  its  solenoidal  (V*b  *  0)  and  lamellar  (V  x  b»  ■  0)  components;  these 
projection  operators  can  be  represented 

jb;  (4 it)-1  V_  x  J  dx  G(x,  x)  x  b  (x,  t) 

1 1  b  !  -  (4n)_1  V  J  dx  G(x,  x)  ?-b  (x,  t) 

and  can  be  used  to  invert  the  differential  equations  x  l>  ■  and  V»l>  «  S 
for  the  field  _b.  For  the  case  of  spacetime  dependent  conductivity,  such 
operations  are  central  to  the  construction  of  the  proper  Hertz  vector,  because 
the  vector  equation  for  A(Z)  one  must  solve  in  order  to  obtain  a  homogeneous 
equation  for  (4b)  is 

V  x  V_  x  A  (Z)  »  c-1  (3  +  4ir  0  )  V  X  V  x  Z  +  V  x  V  x  A  (Z ) .  CB1) 

c 

In  order  to  solve  (Bl)  one  may  define  two  vector  functionals  in  terms 
of  the  Greens  function  G(x,x)  *  G(x-X)  and  the  projection  operators  "• 

A  lamellar  functional  of  any  scalar  is  obtained  from  »,  viz. 

&{s)  =  -  v/|f  G(x-x)  s(lT,t) 

£ 

is  a  curlless  field  solving  V>£  -  s.  A  solenoidal  functional  of  any  vector 
is  obtained  through.!,  i.e., 

2  /f  G(£-Dx AC?.t) 

is  a  divergenceless  field  solving  7  x  V  x  m  X_t.  These  functionals  can 
invert  (Bl)  whenever  the  RHS  is  a  truly  solenoidal  field.  The  orthogonality 
conditions  (11  a,b,c)  and  the  condition  that  7a»VxVxZ  ■  0  are  sufficient  to 
insure  a  solenoidal  RHS  for  (Bl);  one  must  find  the  appropriate  s  and 
For  s  one  possibility  is  the  commutator 


(B  2  a) 


[V*UZ,  "aV-zJ  -  z.Vor  , 

and  for  the  conmutator 

[VxVxcfc,  "5Vx  VxZ_  ]  -  (To*  x  VxZ)  +  V.  x  (Vo-  x  Z)  .  (B2b) 

Taking  the  limit  c  ^V(x,t)  -*■  0  in  (Bl)  one  finds  that  the  field  functional 
A(o)  (Z)  »  c'1C9^  +  4tto)Z  -  /  {—  V^.Z}  -  £C  {— (Vo  x  VxZ)  + 

—  (V  x  73"  x  Z) }  (B3) 

in  a  solution  of  (Bl)  provided  that  the  commutator  (B2b)  is  indeed  a 
solenoidal  field,  a  property  insured  if  and  only  if  Vo'.VxVpcZ  =  VtJiE*  *  0 
as  assumed.  The  result  for  (Z)  gives 

V.A(o)  -  c“1Ot  +  4iro)  V.Z 

Vx7xACo)  -  c“1C3  +  4TO)  VxVxZ 


2  At 

for  all  functions  o(x,t,  Vz(x, t))  because  the  functionals  £  and  £  have 

been  chosen  so  as  to  cancel  all  To"  dependent  terms.  Noting  that  V»A^  has 

now  provided  the  possibility  of  a  proper  gauge  transformation  when  V»£  ^  4> , 

the  remainder  of  (Bl)  must  be  solved  with  a  £fc  functional  if  this  result  is 

to  be  preserved  easily.  Assume  that  A(Z)  ■  A^(Z)  +  £*"  and  the  equation 

determining^  is,  from  (Bl), 


V  *  1  *  & 


HTTP 

c2 


V  x  7  x  [a(o) (Z)  +  £C 


The  complete  solution  is  obtained  by  iterating  for  the  remaining 
unknown  (solenoidal)  component  of  A(Z_) ,  i.e.,  £C.  This  results  in  a 
spatially  ordered  exponential  operator,  i.e.,  (x'.J  >  |x^ |  >  | S3 I • • -  ln 
successive  integrations  (convolutions)  with  G(x,1c) .  The  complete  expansion 
of  A(Z)  is  then  compactly  expressed  as 


B2 


to  all  ordeTS  In  c  V  (x,  t).  With  this  solution  available  one  finds  that 
Z  is  propagated  according  to 

7(7-Z)-7x7xZ  -  c"1  3  A  (Z)  -  0,  CB5) 

—  £ 

which  properly  subsumes  the  relation  (4a)  when  one  identifies  0  •  -  7^  *Z 
and  notes  that  V_  •  A  (Z)  ■  (Z).  In  particular,  (B3)  and  (B4)  imply 

the  novel  gauge  choice 

7  vA  +  c"1  ^3t  +  4ir  o'  (x,  t,  V  (x,  t))]  0  (x,  t)  «  0,  (B6) 

and  the  scalar  potential  evolution  (implied  by  (B5)  or  by  the  substitutions 
for  A(Z),  0(Z)  in  (4a)  )  is  governed  by 

V2  0  -  4*  c’1  Bt  (o(x,  t,  7  (x,  t))  0  (x,  t))  -  c-2  B2  0  -  0.  (B7) 

Prom  these  results  it  can  be  seen  that  the  self-consistent  retention  of 
the  assumed  orthogonality  relations  (11a,  b,  c)  is  possible  when 
o-o  (r,  t)  y  (|  Vr  (r,  t)|  ),  V  -  V*  (r,  t),  and  0  -  -  V  •  Z  -  7-z  C*(*,  t)  *.  0 
This  is  because  the  pure  radial  variation  of  V,  a  allows  the  wave  fields  to 
interact  coherently  with  the  load  over  all  axial  positions  and  because,  cf.(B7), 
the  radial  variation  in  ?  is  Incompatible  with  a  pure  axial  variation  in  <}>  if 
o'  is  non-zero.  The  perfectly  conducting  boundary  condition  on  the  guide  surface 
forces  the  scalar  potential  to  be  attenuated  everywhere  once  a  nonzero  conductiv¬ 
ity  exists  anywhere  in  the  radial  direction.  This  coupling  means  that  only  the 
transverse  waves  (TEM  modes)  arising  through  CC(r,  t)  can  exist  in  the  wave¬ 
guide  formed  by  perfect  conductors.  These  transverse  fields  are  completely 
adequate  for  the  maintenance  of  the  proper  orthogonality,  however,  since 

7  o  *  7  x7  x  A  (Z)  a  0,  7  o  •  7  x7  x  Z  •  0,  V  •  V  -0,  and  7  o  x  V  -  0  are 
—  —  —  —  —  r  —  r 

valid  forever  if  valid  initially. 


r 


The  implication  of  (B3  -  B7)  and  the  self-consistent  retention  of  the 

orthogonality  conditions  is  that  a  proper  Hertz  vector  has  been  chosen.  Again 

the  reduction  to  a  single  propagation  equation  (B5)  for  2  has  simplified  the 

problem,  considerably  extending  the  program  of  reduction  available  in  the 

simpler  case.  In  contrast  to  that  simpler  case  of  arbitrary  but  specified 

p,  J,  the  present  technique  has  absorbed  all  reflected  wave  components  in  the 

% 

non-local  convolutions  involving  G(x,  x)  and  has  eliminated  the  need  for  an 
explicit  source  term  to  produce  them.  In  the  event  that  •  E  cannot  be 
assumed  to  vanish,  e.g.  in  the  case  of  resistive  guide  planes,  this  approach 
can  probably  be  made  to  work  but  the  problem  is  more  difficult  and  remains  under 
study. 

Specializing  completely  to  the  plane  parallel  perfectly-conducting  diode 
geometry  and  using  circular  cylindrical  coordinates,  one  may  integrate  out  the 
trivial  azimuthal  coordinate  in  G(x,7) ,  expand  the  operators  in  (B5)  through 
first  order  in  c  *Vr(r,t)  and  V  Jin?,  and  thus  obtain  the  non-local  (radial) 
propagation  equation  for  Z^r,t)  *  z  CC(r,t).  Suppressing  the  vector  notation 
for  2  and  noting  the  space-time  dependence  of  ?,  this  relation  becomes 


>t  z  -  \  Z  +  c'1  >t  |  Z  +  4ir  c'1/  ir1  Z(rr  t)  »  J  (tj,  t) 

-  4.C-1  to  r  drt  ^  [«tj  Z)  (8^  S)  -  £  Is  [  c'1  J  2  ] 


-  4.TTC 


0» 

1  jt  drl  rl  *"  rl[(V)  *t  *^]9tlZ]| 


Tbfe  observable  fields  Z^  B0  are  then  obtained  from  ZCr,  t),  the  solution  of 
(B8),  as 


-  -  r  *  3r  r  3r  Z  (r.  t) 


(B9a) 


rl 


h  -2. 

c 


and  the  field  (o  E/  x  B.)  then  provides  the  self-consistent  coupling  to  the 
z  o 

fluid  conductor  required  to  calculate  (r,  t). 
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APPENDIX  C 


COUPLING  OP  FLUID  AND  FIELD  EQUATIONS 


If  (B8,  and  B9)  are  to  be  coupled  to  an  annular  fluid  plasma  load, 
represented  by  scalar  fluid  fields  n  ^  Cr,t)  B'Cr  ,t)  ^  T^r ,t)  T^Cr.t) ,  and 
the  convection  vector  field  (r,t) ,  then  one  must  examine  the  detailed 
constitutive  relation  J'  -  oE'  appropriate  to  the  plasma  and  verify  that 
all  plasma  behavior  is  compatible  with  the  orthogonality  constraints  used 
in  the  electrodynamics.  Insofar  as  strict  weakly  coupled  plasma  theory 
exhibits  a  linear  but  anisotropic  Ohm's  law,  the  simple  constitutive  rela¬ 
tion  used  in  developing  the  Hertz  vector  is  not  rigorously  correct.  However, 
it  is  easily  seen  that  to  a  very  good  approximation  the  conductivity  tensor 
is  "diagonal  on  a  slow  timescale",  i.e,  the  radial  E,  J,  can  relax  on  the 
time-scale  characteristic  of  the  electron  e/m  to  cancel  the  radial  currents 
and  to  produce  a  quasi-static  ambipolar  space  charge  and  a  simple  quasi-static 
line-image-charge  radial  electric  field  vector.  In  the  axial  direction  the 
simple  (spacetime  dependent)  o(r,t)  is  also  appropriate  only  on  a  slow  time 
scale  as  the  relative  electron  flow  velocity  field  relaxes  to  a  steady  state 

for  fixed  E' . 

z 

23 

Following  Braginskil  one  may  derive  two  independent  relations  from 
the  kinetic  theory  moment  equation  connecting  d^J  =  y*dsv  ^iVVi^t) 

to  the  impressed  fields  _B  within  the  plasma.  Here  "i"  is  a  species  index 
over  the  appropriate  set  of  one-body  LTE  distribution  functions  for  the 
plasma.  In  the  limit  that  radial  plasma  conduction  currents  are  divergence- 
less  (quasi-static  p  ambipolar),  all  plasma  gradients  are  radial,  and  all 
neutral  convection  is  radially  directed  and  axially  independent,  these 
relations  (with  i-e,  I)  for  the  axial  and  radial  currants  become  raapactivaly , 


"MWH I 


f  •  ■••••*•'  >'•  f- 


f) 


L-X)  ilt  -  -  ^  (1  |e;  -  5,.-1  3rTe  +  ^  J  (Cla) 

V  +  *V  ‘  ^  (<l  +  f^>  Er  +  «"*  "a  x  ®6  +  (i  3rP.  '  3rPl)  "e 


■m  m 

(1  +  2^  -i-su  +  (1  + 
nLj  eTfi  x  z 


®L  e'1  3rT.j 


(Clb) 


where  t 


4  \2t\]  A  e- 


3^  ’ 


U2  "  -  JZ(e  V  * 


x  =  p^Ceg^)-1 


a"i  ‘ 1 " x  * 

and  a,  are  dimensionless  functions  of  ^eTe  as  given  by  Braginskll,  with 
X, x »  denoting  particular  tensor  components.  As  shown  CCla)  represents  an 
equation  of  motion  for  an  "average  electron",  driven  by  and  3rTfi  and 
subjected  to  a  drag  (m  of./ei  )  U  proportional  to  the  velocity  achieved. 

On  a  timescale  characteristic  of  the  electron  inertia  this  "pseudo-particle" 
will  achieve  a  terminal  velocity  determined  by  the  condition  Uz~0. 

Similarly,  (Clb)  represents  a  force  balance  for  an  "average  electron"  which 
is  responding  to  the  space  charge  field  Er  created  by  pressure!  and  temperature 
gradients  and  U  .  Again,  on  an  electron  inertia  timescale  the  apace  charge 
field  self-consistent ly  relaxes  to  a  state  of  equilibrium  stress  in  response 
to  U^»  3rPj»  i a e •  i 


Since  the  inertia  of  the  electron  fluid  is  quite  small  relative  to 


that  of  the  ion  component,  a  fully  self-consistent  treatment  of  the  evolu¬ 
tion  in  Xr(r»c)  which  is  fundamentally  an  ion  fluid  response,  and  Z(r,t) 

(the  electromagnetic  field  incident  on  the  plasma),  is  possible  in  the  approxi¬ 
mation  that  both  U  and  X  relax  rapidly  to  equilibria  values  as  V  ,  Z  evolve. 
For  the  axial  conduction  current  component,  this  approximation  implies  the 
same  constituitive  relation  used  in  developing  Z  in  the  first  place,  with 

a  quasi-static  B„  arising  from  the  thermoelectric  field  -U  n-1  3  T  added 
-  o  x  e  r  e 

to  the  waveguide  field,  of  course.  Substituting  for  J  in  Cla  with  U  -M),  one 

z  z 

has 

*  °E*  -  Jth£  (C2) 

with 

^(r.t)  =  Y  j  J"  £nx(r,t)  Te(r,t)/ax  Cfie(r,t> ,Te(r,t»J  , 


and 


H  m  rw 

Jche  i  6„ 


Me- 


The  conductivity  SXr.t)  defines  the  plasma  location  through  its  density- 

dependent  terms,  if  the  complete  expression  for  T  is  substituted; 

© 

a(r,t)  -|y  T3J2  (r,t)/e2gcTxC£)eTe)A(ne(r,t),  Te(r,t))J; 

(C3) 

and  quite  properly  o  +  0  as  ng  +  0  at  any  radius  r,  because  A  is  driven 
singular. 

Turning  to  the  radial  convection  current  component,  pVr,  the  radial 
field  of  the  linear  image  charge  Er  can  be  expressed  as 


Er<X)  -[/  d?  e  n^.t)  X(?,t)J  •  ^  (C4) 

o 

and  the  self-consistent  x(r>t)  is  then  the  solution  of  the  integral  equa¬ 
tion 


X  h  Vt  ■  r  j(1  +  !f>  [Er(X>  ‘  ^  +  V1  9rTe] 

+  c'1  D  x  B„  +  7-  3  P  -  ^  a  P,  ( 

2  0  n  e  r  e  nL-en  r  A  j 


with  the  fluid  stress  given  by 


(C5) 


fevr.|x  Er(X)  -  Cffl^c)"1  %#Be  -  cahe)_1  3rCPe  +  Pj)  (C6) 

with  m  =  me  +  g-1  nij  and  the  (classically  weak)  viscous  stress  neglected. 

In  a  strict,  completely  self-consistent  electromagnetohydrodynamic 
(EMHD)  field  and  plasma  evolution  in  the  perfectly  conducting  plane  parallel 
waveguide,  one  must  first  solve  (B8)  and  from  (B9) calculate  E  ,  B.; 

Z  w 

then  using  (C3)  determine  U  and  solve  the  integral  equation  (C5)  for  X* 

When  Ez,  B.  and  X  are  available,  the  convection  field  Vr  can  be  evolved  using  (C6) 

and  used  to  correct  the  Z  field  through  (B8) .  To  this  scheme  must  be 

added  a  plasma  energy  equation  for  T  ,  Tj  and  proper  equations  of  state 

for  P_  and  P_. 
e  1 
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APPENDIX  D 


WAVE  EQUATION  QUADRATURES 


If  the  two  dimensional  Green's  function  for  the  wave  equation  is 
applied  to 


3*Z  -  x_13x  x3xZ  -  4ir^f  (E  fjR.Tj) 


(Dl) 


the  formal  solution  can  be  written 

Z(x,T)  -  4 Tfjdtj  dR(2R)|y  dTx  g2(T-T Jr)^  (E^lR.tpj 
o  o  T0 

IT  t-t 


if  /  d<l>/dR(2R)  g2(T-To|R)Z(xoTo) 


(D2) 


o  o 

t r  T-T 


+  h  3T  /  d*/  dR(2R)  82(T-TolR)  Z(Vo) 


o  o 


2u (T-R) 

where  xq  =  xz  +  R2  +  2xR  cos<)>  and  g2(T|R)  -  (Tz_r2)*s  .  Here  u(T-R)  is  the 

unit  step  function  and  T  is  an  initial  time.  This  relation  can  be  the 

o 

basis  for  either  implicit  or  explicit  numerical  methods  because  the 
source  term  involves  the  solution  Z(x,t)  itself.  The  explicit  methods 
derive  from  the  choice  to  expand  (t^)  in  a  Taylor  series  about  Tq  through 
some  fixed  order;  while  the  implicit  methods  arise  upon  noticing  that  ^ 
is  itself  a  time  derivative, 


so  that  a  partial  integration  produces  the 


difference  in  gjlffover  the  time  domain  and  a  convolution  involving  3X  g2> 
The  explicit  methods  are  a  great  deal  less  cumbersjoe  and  more  flexi¬ 
ble  with  regard  to  interactions  with  a  simultaneous  hydrodynamic  calcula¬ 
tion  so  they  will  be  developed  here. 


Dl 


Expanding^  (x  ,  t  +  T)  through  second  order  in  T  i  t-t  (we  suppress 

A  O  O  O 

the  Effnotation  and  use  xQ(x,R,<t))  and  calculating  the  convolution  integral, 
one  finds  ,  J*  |  /  +  (J2  .  R=,>l\ 

Z(x,t)  =  8tt y  d4>y  dR(2R)  jin  ^ - j - J  * 

(4  +  + (£ + r)^c)  -  [n  -  {4  *  I  t+4) 


IT  T 


hf  d*/dR(2R)  TrnF)^  *  J  d*/dR 


(2R)  Z, 


where  the  limit  £  -*■  0  (simultaneously)  in  T+  =  T  +  £,  R+  =  R  +  £  is  under¬ 
stood  and«</  =  3„  *f(x  ,T  ),  etc.  The  limit  arises  from  the  temporal  con- 
T  O  O 

O 

volution  and  step  function  u(T-R)  in  that  one  must  integrate  over  T  just 

beyond  the  6  function  singularity  in  order  to  obtain  a  bounded  contribution. 

The  calculation  of  the  limit  requires  that  one  compute  the  R  integration 

separately  near  R  ->  o,  approximating  ,  ...)  as  constant  over  some 

A  A  0  ° 

small  domain  R  ^  R  <<  T,  Here,  R  may  be  defined  in  practice  as  any  mesh  point 

close  to  the  origin,  such  that  R3  £n«^  <<  1. 

R  o 

2  % 

Performing  the  limit,  making  the  changes  of  variable  R  -►  T(l-v  )  , 
and  x2  •*  a,  and  calculating  the  derivatives  with  respect  to  T  by  Leibniz 
rule  produces  the  final  form  of  the  forward  quadrature  formulae.  For  the 
Hertz  potential  one  finds 


Z<a,To  +  T)  -  16.’|[rH„(|T)  -5f]  '  [e<  +  Ti0  +  ££|  j 
IT  V 

+  871  /d»  fdv  V  In  V  — ■  T2  UQ  +  I«i  +  y-(l  +  fd-v2))^] 
°  77  1 

-  877  fd<t>  JdV  v2T3  [jQ  +  |  ] 

o  0 

7T  1  _ 

+  —  Jd<b  fdv  [zo  +  TZq  +  TZ^*  (2  cos<J>  Va(l-v2)  +  2T(l-v2))] 


and  for  its  partial  time  derivative, 

D2 


( 


Z(a,To  +  1)  -  16ir=j[ft=  ta(£)  -  +  T_4]  + 

[n  +®r]-t<+  t4  +  ItH]|  + 

A 

7f  V  _ 

Sir /*/ dv  v  fcn  ^2T^jf  +  3T2^q  +  2T3(1  +  A  (l-v2))^)^ 


2(T2cos4>  /a(l-v2)  +  T3(l-v2)) -Ce>£'  +  T^/'  +  T2(l  +  4(l-v2))i 

O  o  2 

IT  1 

-  8tt  f d<t>  f dv  v2  3T2*/  +  3T3,ef  +  (2jgl'  +  hd')  (T^l-v2) 

J  J  (  o  o  o  2  o 


o  o 


+  7a(l-v2)  T 3  cos4>)  | 


^  r' 

'  dv 


Z  +  2T(l-v2)Z '  +  2(2Z '  +  TZ')  (cos<J>  V a(l-v2) 
o  o  o  o  v 


+  i?  /d<*>  /« 

o  o 

+  T(1-v2)^+  4  (cos  a(l-v2)  +  T(l-v2) ) 2  TZ^ ' 


with  the  notation 


z'  =3  Z(a  ,  t  ),  Z’  =  3  3  Z(a  ,  t  ) 

o  a  o  o  o  at  o  *  Cr 


<  5  3a.4(Vo)  *  4'  E  8  a  \  ^o  (ao*To) 


-  -  o  a 

o  0  0 


H  8a.3^o(ao  -V* 


and  the  arguments 


a  *  a  +  2T  cos<J>  /a(l-v2)  +  T  (1-v2) 


The  quadratures  can  be  carried  to  higher  order  partial  time  deriva¬ 
tives,  but  these  two  constitute  the  minimum  set  required  to  specify  the 
source  term  4*£(a,To  +  T)  for  a  subsequent  time  step.  The  quadratures 

over  0  and  v  can  be  represented  accurately  using  well-known  integration 
24  25 

formulae  ’  over  the  unit  circle  once  the  discrete  representations  for 
ZCa^T^)  and  Z(ajT^)  are  interpolated. 
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